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Introduction 

The work performed under NASA grant #N3R 14-004-008 
(originally NSG-694) and its extensions cohered many aspects 
of the fluid dynamics of coaxial flows and of rotating flows. 
Previous results of this work received publication in the form 
of both NASA contractor reports and professional journal 
articles. Lists of these reports and papers are presented in 
an Appendix of this report. In the body of this report two 
different analytical investigations of fluid dynamics problems 
of relevance to the gaseous core nuclear reactor program are 
presented. Since all prior work had already been published, 
the results of these investigations constitute the final report 
of work done under the Grant. 

The first of these two works, "Numerical Solutions of 
Driven Vortices of Binary Fluid in Cylindrical Geometry" is 
an analysis of the vortex type flow which appears in the nuclear 
"light bulb" concept. In this work, a numerical treatment is 
developed for the rotating flow which includes a description 
of the nuclear fuel addition. The problem is formulated with 
the complete Navier-Stokes equations in order to show the 
interaction between the fuel addition, the main flow and the 
boundary layer flow in an accurate manner. The results pre- 
sented are the first to show holdup of the nuclear fuel for 
the case of steady fuel addition. The results for fuel holdup 
are discouraging and it remains to be seen whether optimisation 
of injection location will provide substantial improvement. 

The second work, "Laminar Confined Coaxial Entrance Flow 


r 


4 


with Heat Generation/' is an analysis of the fluid flow in 
the fuel inlet region for the coaxial flow gaseous core nuclear 
reactor concept. This is a region in which analysis is 
difficult due to the existence of very large gradients. In 
the analysis , the boundary layer form of the equations of 
motion is used because the large gradients make a more exact 
formulation extremely difficult to solve numerically. Results 
demonstrate that the large expansion due to temperature change 
will take place axially unless strong radial pressure gradients 
are developed by induced non-axial velocity components. This 
represents a significant difference from the present conceptual 

description. 

These two investigations are the culmination of the NASA 
sponsored I IT work on gaseous core nuclear reactors. The 
first of the two investigations extends the development of 
numerical methods for the solution of the Navier-Stokes 
equations for appropriate geometries to the case of rotating 
flows and almost completes the gas core program requirements 
in this area. It is still necessary to extend the capability 
for Navier-Stokes equations integration to spherical coordinates 
and to increase the range of Reynolds numbers for which 
solutions can be found with reasonable computer time usage. 

The second investigation demonstrates that even the 
conceptual design of the coaxial flow reactor needs further 
development. The extension of analytical capability to spherical 
coordinates will greatly aid in the further conceptual develop- 
ment of the reactor. 
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ABSTRACT 


’Ph® steady state laminar notion of a viscous , incom— 
pressible and binary fluid is studied for a rotating flow 
in a cylindrical geometry. The mathematical model employed 
is a cold flow simulation of the fluid mechanics of the 
Light-Bulb concept of the gaseous core nuclear engine. The 
radial inward convection of angular momentum of the coolant 
gas through the porous periphery drives the vortex, mixes 
with the fuel gas and eventually leaves the chamber as a 
mixture through a small exit hole on one end wall. The fuel 
gas addition is mathematically simulated by an interior mass 
source . 

A new computation method is presented to account for 
the interior mass addition. The method uses a po te ntial 
function in addition to the stream function and vorticity 
used in current methods for solving viscous incompressible 
flows. By assuming that the fluid's velocity consists of 
two parts, one that is the gradient of a potential function, 
and the other calculated by appropriate derivatives of the 
stream function, it was found possible to satisfy the con- 
tinuity equation which contains a source term. This method 
can also be applied to unsteady compressible flows where 
wis unsteady continuity term is viewed as a source term 
while iterating over an intermediate time level. This 
technique is applied to driven vortex flows of a binary 
system and carried through to numerical results. 
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The governing equations are formulated in terms of 
Stokes' stream function , the tangential vorticity and a 
potential function as a time dependent problem. Steady 
state solutions are then obtained as large time solutions 
of the time dependent problem. A numerical method by finite 
differences , specifically the alternating-direction-implicit 
method, is employed for solving the equations. Once the 
velocity field is solved, the pressure is obtained by 
solving the Poisson equation of pressure with Neumann 
boundary conditions that are determined by given kinematics. 
Por purposes of computational efficiency, the solution 
method by Fourier analysis is based on fast Fourier trans- 
forms. 

Numerical results were obtained for axial Reynolds 
numbers of 1 and 20, tangential Reynolds numbers up to 100 
and aspect ratios of 1 and 2. Without rotation, the primary 
flow was found to be more or less uniform through the bulk 
of the chamber. However, with rotation the bulk of the 
fluid flows through the relatively narrow region of the 
boundary layers formed along the stationary end walls. Also, 
the flow in the main body of the chamber consists of 
recirculating secondary flow cells. The intensity of this 
flow feature increases with increasing rotation for fixed 
through-put and decreases with increasing through-put for 
fixed rotation. The results obtained for binary flows show 
that the hold-up of internally injected mass increases with 
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increasing rotation for fixed primary flow and Schmidt 
number, finally , the containment factor increases with 
decreasing the density ratio ( density of fisrionable mat- 
erials to density of coolant ) with other parameters kept 
constant. 
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NOMENCLATURE 




In the following, symbols and abbreviations used in this 
thesis are listed. The symbols used only locally are ex- 
plained there and hence not included here. 


English letters: 

a aspect ratio 

b non-dimensional geometry parameter as defined 

in (59) 


AB 

D/Dt 

e (r,z) 



A-B 


h 


H 

k 

L 


M 


n. 


n. 


B 


N 


P 


mass fraction-of— component A 

mass diffusivity in binary diffusion 

material derivative 

local error function in the continuity 
equation 

discretization error 

fractional mass addition of component A 
relative to that of component B 
Ar 

vorticity-like function defined in Eq. 15 
Az 

length of vortex chamber 

number of mesh spacings in z coordinate 

mass flux vector of component A 

mass flux vector of component B 

number of mesh spacings in r coordinate 

static pressure 
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<i. 


‘s 


R 

Re 

Re^ 

Sc 

t 

u 


u 

V 

-► 

V 


inj 


w 


ex 


m 


z 

Greek Letters t 
a 
@ 
r 

r_ 


iteration parameter in elliptic ADI equation 
of potential function 

iteration parameter in elliptic ADI equation 
of stream function 
radial coordinate 

interior mass source of component A 

exhaust hole radius 

chamber radius 

axial Reynolds number 

tangential Reynolds number 

Schmidt number 

time coordinate 

radial component velocity 

radial injection velocity on periphery 

tangential component velocity 

velocity vector 

tangential velocity on periphery 

axial component velocity 

axial velocity at the exhaust hole 

reference velocity for flow due to interior 

mass addition as defined in (60) 

axial coordinate 

At/hk 

hA 

angular momentum 

angular momentum on periphery 
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A 

Ar,At,Az 


V 

V 2 

V 
Vx 
e 
C 
n 
x 


V 


p 



p 


B,P 


4> 

* 


$ 


(it 

Subscripts : 


A 

B 


6 

1 

j 


hk 

radial mesh size, time step, axial mesh size 

respectively j 

del operator 

Laplacian operator 

divergence operator 

curl operator 

, 

potential function parameter defined in (59) 
tangential component vorticity 
containment factor 
volumetric coefficient of viscosity 
shear coefficient of viscosity 
mixture density 

density of component A in mixture 

i 

density of pure component A 
density of component B in mixture 

i 

density of pure component B j 

potential function j 

■1 

Stokes' stream function ^ 

| 

vector potential 

I 

vorticity vector j 

refers to component A 1 

j 

refers to component B 

refers to exhaust hole 

mesh point in radial direction 

j 

mesh point in axial direction s 
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Superscripts : 

n refers to time level 

* refers to non-dimensional variable 

Abbreviations : 

ADI alternating direction implicit 

FFT fast Fourier transform 

SOR successive over relaxation 
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CHAPTER I 
INTRODUCTION 


Rotating flows have long been a subject of investiga- 
tion to fluid dynamicists as applied to Rangue-Hilsch 
tubes cyclone separators , heat exchangers , MHD vortex 
power generators and many other devices. Recently ad- 
vanced concepts of gaseous core nuclear reactors have fur- 
ther stimulated the study of vortex flows. The subject has 

( 14 ) 

also been studied for purely academic purposes. 

The present study has directly been motivated by one 
of the advanced reactor concepts, i.e. the Light-Bulb con- 
cept , and is concerned with containment of fissionable ma- 
terial and how the phenomena involved in vortex flows in- 
fluence containment. 

Definition of the Problem 

In the Light-Bulb reactor the fissionable material 
is kept inside a cylindrical container and the propellant 
flows outside of the cylinder tnrough an annular space. A 
transparent cylindrical shell separates the propellant from 
the fissionable material. The fuel is maintained in a vapor 
state at very high temperature, and transfers heat to the 
propellant through the transparent separating shell by 
radiation. To maintain the integrity of the shell which 
partitions the fissionable material and the propellant, it 


* For all numbered references, see bibliography. 
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is necessary to inject coolant gas at the shell. *?he coolant 
gas is intended to play another important role. Since 
the coolant is injected with rotation , a vortex is formed 
which then causes the formation of boundary layers on the 
end walls of the chamber. The flow of fluid in and out of 
the boundary layers induces recirculating secondary flow 
cells in the main body and fluid dynamically retains the 
fissionable material in this region of closed flow cells. 

The coolant gas fed to the vortex chamber is withdrawn through 
an end wall, cooled and then recirculated. Fuel entrained 
by the coolant is separated before recirculation. 

In actual operations, a required amount of fissionable 
material called the critical mass has to be maintained in 
the reactor chamber. This necessitates an attainment of a 
certain fuel concentration. The coolant gas that is in- 
jected to cool the wall separating the propellant gas and 
fissionable material entrains some of the fissionable materi- 
al while passing through the chamber. Although the required 

concentration level of the fissionable material can-be- 

maintained by replacing the amount entrained by the coolant 
gas it is very desirable to keep this as small as possible. 

By forming a vortex and therefore inducing recirculating 
flows, the critical fuel mass may be held in the closed flow 
region with a minimum loss. Because of the particular fea- 
tures of the confined vortex flow, the manner in which the 
fissionable material is introduced is important. It has 
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been suggested that the fissionable material be injected 
from an end wall in the form of small pellets so that it can 
reach the retention region without being swept out of the 
chamber by the boundary layer flows. 

The present study consists of a numerical solution of 
a binary fluid system in a cylindrical container of finite 
length as shown in Fig. 1. One Of the components can be the 
fissionable material and the other coolant. Addition of 
fissionable material is simulated by an interior mass source 
in the annular containment region. The coolant fluid enters 
through a rotating porous periphery and leaves axially 
through a small exhaust hole on one end wall. 

Because of the particular manner in which the fission- 
able material is introduced, the continuity equation con- 
tains a source term. For this reason the conventional com- 
putational scheme for obtaining numerical solutions in terms 
of the stream function and the vorticity fails. 

A new computational method is developed here to treat 
this problem of interior mass addition by means of a con- 
venient formulation which consists of the stream function, 
the vorticity and one other function. This function is a 
potential function which is introduced into the formulation 
so that the continuity equation that contains a source term 
is satisfied. It may be mentioned that the same technique 
can be applied to time dependent compressible flow problems. 
In such problems the unsteady term in the continuity 


uniform radial injection of simulated coolant 
porous rotating periphery 

•simulated fuel addition 
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exhaust 

hole 


exit flow 


stationary solid walls on both sides 


Pig. 1. The Plow Configuration and the Adopted 
Coordinate System. 



equation Is treated as a source term in the iteration pro- 
cess over the time step. The method is used in this study 

to obtain numerical solutions for binary driven vortex flows. 

The basic differential equations are written in un- 
steady form in terms of the stream function, the vorticity 
and the newly introduced potential function. A numerical 
method by finite differences has been employed as the solu- 
tion method. The altemating-direction-implicit met’.'od is 
used for both the parabolic and elliptic equations. After 
the kinematics are solved, the pressure is Obtained by solv- 
ing the Poisson equation of pressure. The method of Fourier 
analysis is used for this purpose in which the fast Fourier 
transforms have been employed to obtain an efficient comput- 
ational procedure. 

Literature Review 

Rather than attempt to give an elaborate historical 
survey of vortex flow investigations the discussion here is 
limited only to those that are directly related to the pre- 
sent problem. Earlier analytical approaches to vortex flow 
problems assumed the flow to be one-dimensional J 11 * There 
have been other investigations which considered the boundary 
layer flow formed over a finite disk with some specified 
outer velocity distribution. These studies have 
not considered the influence of one part of the flow on the 
other, i.e. the bulk flow affecting the boundary layer flow 
and vice versa. Accordingly, results of these studies do 
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not render sufficient information to describe the confined 

vortex flow problem. The importance of the interaction be- 

( 2 ) 

tween the two flow regions was first noted by Anderson; 
who considered the laminar boundary layer over the station- 
ary end wall. More recently # Rosenzweig # Lewellen and 
Ross* 31 * studied the interaction problem analytically# in 
which the flow was divided into three subregions# the end- 
wall boundary layer# the primary flow and the central vortex 
core. For high rotation rates# they obtained a similarity 
solution for an assumed turbulent boundary layer. A com- 
prehensive review on confined vortex flows can be found in 
Lewellen* 22 * ' * 24 * and Rott and Lewellen.* 33 * 

Experimental work on confined vortex flows includes 
the flow visualization studies of Roschke* 3 *** and Rossi 32 * 

In both of these the secondary flow pattern shown in the 
visualization results was obtained for jet-driven vortices. 
Roschke studied also the effects of aspect ratio and the end' 
wall geometry on resulting flow configuration. 

However# no numerical work can be found in the litera- 
ture that considers the driven vortex flows in the config- 
uration which simulates the Light-Bulb reactor concept. A 
numerical method applied to the flow configuration which is 
probably closest to that of the present study is given by 
Hombeck. * 17 * He only considered the limiting case of 
small aspect ratio in which the recirculating flow could 
not occur. Nevertheless# the literature of rotating flows 







A22 


in a system completely enclosed by soxid boundaries or in 
a semi-infinite or infinite region is large. Equations of 
hydrodynamics may be written in terms of velocity and pres- 
sure or in stream function and vorticity. Of the many 
studies that employed the stream function-vorticity form- 
ulation, the important ones include the work of PaoJ 25 * 
Briley Lava nj 2 ^ and Aziz.^ Pao considered the second- 
ary flows in a finite length cylinder with solid boundaries. 
One end wall was stationary and the other wall rotated with 
the periphery. Briley investigated the time dependent 
motion of spin-up, spin-down and revers e- rotation s . T.avan 
examined the development of swirling flows in a circular 
duct. On the other hand, in a few cases the velocity-pres- 
sure formulation has also been used for numerical solutions 
of rotating flows. Williams considered a three-dimen- 
sional thermal convection problem with rotation in an an- 
nular region. It is the author's belief that the work of 
Williams is the only published one that uses successfully 
the velocity-pressure formulation in rotating flows. Final- 
ly, it may be mentioned that no publications can be found 
that consider the rotating flow of a binary fluid system. 

It is the belief of the author of this thesis that the 
method developed here is the one most suitable for obtain- 
ing numerical results for confined vortex flow problems 
of binary fluids. 


CHAPTER II 

MATHEMATICAL FORMULATION 

In this chapter a set of differential equations is 
written for a binary fluid system. One species of the 
system is introduced into the chamber as an interior mass 
source, whereas the second component enters through the 
porous periphery. The equations are obtained from the 
conservation laws of species, momentum, energy and the con- 
stitutive laws. The result thus obtained yields a set of 
equations in terms of velocity and pressure. By eliminating 
pressure from above velocity-pressure formulation, a 
second formulation is derived in which the equations are 
expressed in terms of stream function, vorticity and 
potential function. The merit of the stream function- 
vorticity formulation in numerical integration schemes is 
discussed in this chapter. 

Underlying Assumptions 

The important assumptions made in this study are 
listed below. The basic differential equations are de- 
rived under these conditions. 

(1) Steady, isothermal, incompressible and axi-sym- 
metric flow. 

(2) Binary fluid system. 

(3) The first species enters the chamber as an 
interior mass source with negligible momentum. 
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(4) Density varies with composition only. 

(5) Constant viscosity. 

(6) Mass diffusion is assumed to take place solely 
by concentration gradients. 

The compressibility effect is neglected because the vari- 
ation of pressure in the system is small when compared to 
the high absolute pressure. Constant viscosity is assumed 
because the objective of this study is to model a cold flow 
simulation of the Light Bulb reactor. And the effects of 
concentration variations on the viscosity are small. The 
viscosity would of course not be constant in the reactor 
itself because of the large temperature variations that 
exist there. Other assumptions that are considered minor 
will be discussed in the specific places where they occur. 
Governing Equations 

Diffusion and Continuity . The continuity equations 

for components A and B are obtained respectively by taking 

(7) 

the mass balances of individual components: 


8p & 


+ V»n, 


( 1 ) 


^ *► 

7T * v ‘ n B " 0 


(2) 


Adding these two equations yields the continuity equation 
of the mixture: 




O' 1 ' o 


* >< l*ot T\\l 

-.. . : " ’ /(*- “■. 
« „ •! 









|£ + ?• (pv) = r A 


(3) 


The mass diffusion equation is derived taking into account 
the effect of the concentration gradient only. By sub- 
stituting Fick’s law of binary diffusion into the equation 
for species A, Eq. 1, and making use of the mixture con- 
tinuity, Eq. 3, the species A equation can be expressed 
as follows when the mass diffusivity is constant: 


DC, 


p Sr ’ D AB 7- tp7o A> + r A (1 - C A> 


(4) 


Momentum Equation . The momentum equation with con- 
stant viscosity in the absence of body forces and for 


(3) 

Newtonian fluids is: 


g| = -Vp +yV 2 v + (Up)V(V-v) 


(5) 


Dt 
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With Stokes' hypothesis, 3X + 2y » 0, the momentum equation 
reduces to: 


pjjv = „7 p + y y 2 v + juV (V*v) 


( 6 ) 


'Dt 


The tangential component of this vector equation for an 
axi-symmetric flow is: 






$ 


1/,.I1- (if.M tH,u.|i,n.iK.. - ;1U»M< ♦ , *'H> tw * ( ¥ «' >■ l-»« r l 











[K*f] -“[^-7] 


(7) 


where the operator under axial symmetry is: 


§t ' W + u Sr + 4s 


This tangential component velocity equation is expressed 
in terms of the angular momentum by the variable change 
r * rv: 


«-*[ 


r - — r + r 

1 rr r r zz 


]• 


( 8 ) 
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The equations given above constitute the velocity- 
pressure formulation and numerical solutions for the veloc- 
ity components, density and pressure can be obtained from 
their finite difference analogue without further modifi- 
cation. However, difficulties have been reported in the 
literature in solving the hydrodynamic equations when 
the pressure is included. Therefore in the formulation 
used here, the pressure is first eliminated by introducing 
the vorticity and the vorticity is then expressed in terms 
of the stream function. In this manner the same problem 


* The subscripts in the equation represent differentiations. 
This convention is extensively used throughout this 
thesis. 
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is reformulated in terms of the stream function and vor- 
ticity. Since this formulation cam be solved without much 
difficulties it is adopted in this study and will be de- 
veloped in the following. 

Vorticity Transport Equation . The vorticity trans- 
port equation is obtained by taking a curl of the momentum 
equation which for constant viscosity and no body forces 
is expressed as: 

p [§1E + “ (w* V) v j + Vp x ~ = iiV 2 u> ( 9 ) 

For axi-symmetric flow it is sufficient to consider only 
the tangential component of the vorticity vector: 


’[bI* <“ r + w 2 K - i(v\] ♦P z <5|-£> - 


P ~ 
p r Dt 


y(v 2 c - £*) 


7 ' 


( 10 ) 


A new function is introduced next by a change of variables: 


H * 


£ 

r 


( 11 ) 


This changes the anti -symmetric function 5 to a symmetric 
function H with respect to r, and upon substitution Eq. 10 
reduces to a simpler equation that is especially convenient 
for numerical integration. 
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p [§!* (v-5)h - ij ( v 2 ) z ] + ^<§|-£> 

M(7 2 H + |B r ) (12). 

An important simplification that results from this change 
of variables is that the term, C/r^ no longer appears in 
the equation. 

Potential Function Equation . A potential function 
of an auxiliary nature is introduced in such a way that 
it identically satisfies the continuity equation which 
has a source term. In more conventional problems of in- 
compressible flows or compressible steady flows, the Stokes' 
stream function satisfies the continuity equation which 
in that case has no source term. 

In order to make the discussion general, a vector 
potential is used, which for an axi-symmetric flow reduces 
to a single non-vanishing component, i.e. the tangential 
component . Let 

Pv = V X | - pV$ (13) 

With this definition of the mass average velocity, the 
continuity equation for steady flow becomes: 


* 


Eq. 12 is called the vorticity equation hereafter. 
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7* (7 x $ - p7$) = r A 
Since 7» (7 x $) =0: 

7* (p7$) = -r A (14) 

This is the equation that the potential function must 
satisfy which together with the vector potential J yields 
the velocity vector through Eq. 13. When an unsteady flow 
of variable density is considered, the term 3 p/3 1 can be 
treated as part of the source term r A . 

The potential function Eq. 14 and the expression for the 
velocity vector Eq. 13 are expanded as in the following 
for axi-symmetric flow: 


i p y P„ r. 

a + (£ + _£ )d> + a + + _£1 

y rr r p ' y r y zz pz p 


(15) 


and 


u 


" pF*z " 4> r ; w = £r*r “ ♦ 


(16) 


Stream Function Equation . The defining equation 


of vorticity is: 





Upon substitution for v from Eq. 13 and using the vector 
identity, V x (V$) = 0, the vorticity is expressed as: 

w = 7 x (^7 x J) (17) 

The tangential component of Eq. 17 is after change of 
variables: 


(i_ w, ) + iu ) + 

l pr y r'x l pr y z'z 


rH 


(18) 


Mixture Density . The mixture density neglecting 
effects of pressure changes is a function of composition 
(mass fraction) only: 


P 



P 


Ail 


+ c A (l 


P B,P 


(19) 


Pressure Equation . Two different solution methods 
are considered for the pressure. First, the pressure can 
be obtained by integrating the momentum equations. A 
second method is to solve the Poisson equation obtained by 
taking the divergence of the momentum equation. Both of 
these methods are used in the work reported here. 

In the stream function-vorticity formulation which 
is used throughout this work, the kinematics are solved 
first independently of pressure. Therefore, the pressure 


A31 






may be calculated from known kinematics by either method. 

The momentum equation is rewritten in a form con- 
venient for the line integration of pressure: 

* “P§£ + + ^y7(7»v) (20) 

Also the Poisson equation of pressure obtained as indicated 
above is: 

V 2 p « 7* + y7 2 v + |y7(7«v)J (21) 

It may be mentioned that when the governing equations are 
formulated in terms of the stream function and vorticity , 
the error in pressure calculation does not affect the 
solution of the kinematics. Primarily this is the reason 
why the stream function-vorticity formulation has been 
more extensively used than the velocity-pressure formul- 
ation. Aziz ^ ^ in his studies with both formulations found 
the one using the stream function and vorticity advantag- 
eous even for three dimensional problems in which he had 
to calculate all three components of both the vorticity 
and vector potential (stream function) . 

Transpo rt Equations in a Conservation Form . In 
obtaining numerical solutions by finite differences, it is 
preferable to express the convective terms in a conserv- 
ation form. The basic idea# as can be found in Emmons 


M 


*.. ■ 


1 
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is to conserve the transport quantity even with finite dif- 
ferences. The convective terms expressed in terms of a 
general transport quantity , q, have the form: 

(pv*V)q * V* (pvq) - q(V-pv) 

Using continuity, Eq. 3, under steady .-tate conditions the 
above equation can be expressed in the following con- 
servative form: 

(pv*V)q = V* (pvq) - r A q (22) 

This equation in cylindrical coordinates under axi-sym- 
metry reads: 

p(uq r + wq 2 ) * |(rpuq) r + (pwq) z - r ft q (23) 

Boundary Conditions 

Since all of the differential equations are elliptic 
in the spatial variables, boundary conditions on the 
entire boundary are required. Along the center line the 
condition of symmetry is prescribed, whereas on solid 
boundaries the no slip condition on velocities and a van- 
ishing concentration gradient are imposed. On the porous 
peripheral wall, a uniform injection velocity and zero 
flux of species A are assumed. And at the exhaust hole a 
porous plug condition and zero concentration gradient are 


specified. Under these conditions, the total mass flow 
rate through the exhaust hole is not known a priori. 
Therefore, at each time step of the numerical computation 
the total mass flow is computed from the available solution, 
and from the overall mass balance the injection velocity 
at the periphery is determined for the boundary conditions 
of the • time step solution. 

Th. tundary conditions for the stream function are: 


tj> (0 , z) ■ 0 (24) 

tHr,0) * 0 (25) 

On the periphery, u * ” u inj = constant, and hence from 
Eq. 16: 

*z m ■° R( * u inj + V 


And upon integration: 

z 

i|>(R,z) = -R ^ p *" u inj * 


(26) 


Also 


<Hr,L) - t|<(R,L) for r > r e 


( 27 ) 


At the exhaust hole, from Eq. 16: 


.... f , ' • '^: 1 • 77^ •►••• •• • •v^-^yi < ~v*r».*fe»» .y' 
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i|» (r,L) - / Pr(w ex + 4> z )dz for r < r 


e 


( 28 ) 


The boundary conditions for vorticity are obtained by 
evaluating the stream function equation on the boundary for 
no slip conditions. Especially with the introduction of 
the potential function which by nature allows slip 
velocities, the no slip conditions are attained through 
the boundary conditions on the vorticity. 

With respect to the axis, r * 0, the axial velocity 
is symmetric whereas the radial velocity is anti— symmetric , 
and hence this vorticity like function H is- symmetric. 


H r (0,z) = 0 


(29) 


At the solid boundary, z ** 0, w *-0 and u = 0. Hence 


H(r,0) 


- I(i-d) ) 

r'prz z 


(30a) 


in which the no slip condition, u * 0 is satisfied by: 


<J > 2 - -Pr<fr r 


On the periphery, r = R, u = -u in j, w » 0 and ♦ * 0. 


(30b) 


H(R,z) 


- !(!-* ) 
r prr r 


Hence 


- r . • 


! T •• ,n ■ 


M. iwp w «Wf iB iai 


’ o v, Q 


* ,.. 5 *' *, » . M .. ..* ’« • m • \ n 


M i tM | tlfltll ijk i £« ■ 1 1 


Expanding this and applying the no slip cor.di lon, w - 0 
results in: 


H(R,z) ® - -^r '!> r 
pr 


Ola) 


where the no slip condition, w = 0 is ensured by an add- 
itional condition: 

tli =0 

y r 


(31b) 
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At z * L, u * 0 and hence 


H(r,U - - ;[<jhVr + ‘pr^'t] 

subject to the condition t|> 2 * -pr<j» r . 
Expanding the above expression gives: 


H(r,L) = - jrj^^r + pr ( " p Z ^ z + * zz 

which is subject to the following condition 
u * 0: 


)J (32a) 

to satisfy 




-prij), 


(32b) 


Boundary conditions for the angular momentum are 


r( 0 ,z) « 0 
r(r,0) = 0 
r (R,z) = r o 


(33) 

(34) 

(35) 


i 

3 

j 

| 

i 

j 

i 

i 

i 

i 

i 
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T(r,L) * 0 (36) 

No mass flux of component A across the periphery and 
zero concentration gradient at the exhaust hole are assumed. 
Along the center line the condition of symmetry gives:* 

c (0,z) * 0 (37) 

T 

Across solid boundaries the total mass transfer as well as 
the velocity vanishes. Hence from the expression for the 
mass flux, 

n A * pvc - pD^Vc 

the following boundary conditions result: 

c_(r,0) = 0 (38) 

2 

c z (r,L) * 0 (39) 

At the periphery the condition of no flux of component A 
yields : 

D AB c r (R r z) + u i n j°( R ' 2 ) “ 0 

The boundary conditions on component B can be deduced from 
those on component A. The mass fluxes of each component 
are: 

* The subscript A is suppressed hereafter in the mass 
fraction of component A. 
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n A * P^ c b " °AB an” 


n B - pvc A - d ab an” 


where n represents a local coordinate normal to the boundary 
Adding both equations yields: 


“a + ^ * * 


At the periphery, where n A = 0s 


*B- * \ 

i 

i 

i 

) 

The boundary conditions for the potential function 
are the symmetric condition on the axis and zero normal . 

gradient on solid boundaries. However, at the periphery j 

and exhaust hole an- arbitrary functional value, « - 0 is I 

assigned. In this manner, the conditions on velocity are I 

satisfied by assigning appropriate values to the stream 
function on these boundaries. 


<J> r (0,z) * 0 


♦ z (r,0) * 0 




A3 8 


<t>(R,z) <s 0 

(43) 

$ z (r,L) = 0 for < r < R 

(44a) 

4>(r,L) = 0 for 0 < r < r _ 

« « © 

(44b) 

dimensional Formulation 



In this section the system of equations and boundary 
conditions in the stream f unction- vorticity formulation as 
obtained in previous sections are non-dimensionalized with 
properly chosen reference quantities. The characteristic 
reference values are chosen such that the normalized vari- 
ables remain of order one. This is important since only 
then are the non-dimensional governing parameters meaning- 
ful as they appear in the non-dimensionalized governing 
( 22 ) 

differential equations. ' The non-dimensional variables 
are written with a superscript (*) , whereas the reference 
quantities are written with a subscript (ref) unless al- 
ready specified. Let 

r* * r/R, z* * z/L, t* = t/t f 
u* = u/w , v* * v/v , w* = w/w 

P* = P/P B(P . ** - ♦/♦ref' H * * "/“ref (45> 

r* * r / r o ' ** ■ */*ref' r *A ■ V r A,ref 


p* * p/p 


(45) 


ref 

Substitute the above relations into Eq. 16 and obtain: 

b 2 

U * 2ap*r* “ e< ^r* (46a) 

W * ° 2p*r* ♦£* " §♦%* (46b) 

In above equations two reference quantities are determined 
and three non-dimensional parameters are introduced: 

*ref = 7 p B,P r e w ex ; ♦ref * w m R 

b = r e /R; a = L/R? e = w m /w ex 

where w m represents the reference velocity for interior 
mass addition as the total mass addition averaged over the 
entire periphery. 

Substituting (45) into the stream function Eq. 18: 

(pip- + Vpip " " V* H * ( 47 > 

a b 

“ here H ref ■ *ref / ‘ > B,P r e ia U8e<J - 

The potential function Eq. 15 is normalized using (45): 


5 


I 1 
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^r*r* + ( F* + p*“ )<(, J* + 75^5*2* + 7? pT^J* + ^ “ 0 


a 


(48) 


where r. = p_ _a ,/r 2 . 

A, ref K B,P y ref' 

Now the vorticity Eq. 12 is normalized to: 


H t* * ?f[ H J*r* + r**r» + Vl* 2 * * gfr<r*p»u*H*) 

V» ft 


- K<p* w * b *>z*] + E 2 [ £ pr - <»*. + j£ + iw*,)J 


H* 


2b 


Rer 1 


fRe| 2 Re 2 p** , 

*** ' + s"*“S 


z* 


Re 

Re 


t r» 2 . Re2 p£* i 1 

7 + — p T <u»w»» + £**w**> J 


(49) 


where Re = P B , p w ex r e / M , Re t - P BfP v 0 r e /p and 
fc ref = p b,P r2/m * 

The angular momentum Eq. 8 after non-dimensionalization 
reduces to: 

r t* “ p T [ r ^* r * ~ F^r* + “7 r z*z* ” §§* (r*p*u*r*) r * 

^ A 


" §§(p # w*r*) 


-] 


ReerJ 

bp* 


•r* 


(50) 
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The equation for diffusion, Eq. 4, is normalized to: 


't* 


if . ,1 . p r*. . 1 p z» "1 

g; [ c r*r* + ( ?f + p*" ,c r* + ^ P* -0 **] 


Reer£ 


“ EfSr*< r * (> * u * c * ) r* " 5^f (p * w * c) z* + ”Ep*~ (51) 


where S 


c “ P B,P°AB * 


The mixture density, Eq. 19, is non-dimensionalized toj 


5* a 


P A,P /p B,P 


p A,p/ p B,P + ° 


“ P A,P^ P B,P^ 


(52) 


Finally, Eq. 20 and 21 involved in pressure calculations 
are non-dimensionalized to: 


Pf* “ 


-P* (u«uf + £ uj. - $r-) + 


Re t v* 2 , . b 


u* 


+ £e ixi b + £ + bz^T* 


(53a) 


1 * 

a?z* 


-p* (u*w£* + £w**> + |^V* 2 w* + 


+ £ + k*>z* 


(53b) 



where the reference pressure is determined: 
p ref * P B,P w ex 

Non-dimensional governing parameters and reference 
quantities determined above are summarized below: 
Non-dimensional governing parameters: 

Geometry parameters; 

a = L/R, b = r fi /R 
Potential function parameter# 


and 


e ■ w n/ w ex 


Re = p B,P w e xV 1 


Re t * °B,P v o r e /u 


( 55 ) 


S 


C P B#P D AB 


Reference quantities: 


‘ref * p b,p r2 /''- *ref * I’b.p'Kx' 


I 


*■* , .. 


i -t 1 


AliiUWkl 
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mu 
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H ref * ^ref^ p B,P r e # 

( 56 ) 

*ref " w m R ' r A,ref " p B,P*ref /R2 ' 
w m “ f A-B u in j' r o * Rv o 

The boundary conditions are non-dimensionalized in the 
same way and presented in the summary section of the dif- 
ferential equations and the boundary conditions in 
Chapter III. 


wmrrA wjMJM-ff hi m-w. vmt *m vm um hjmpusjjfij iwm$ 
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CHAPTER III 
METHOD OP SOLUTION 

Due to the extreme complexity and added difficulties 
arising from the non-linearities of the equations, any hope 
for an exact analytical solution is remote. Therefore, it 
was decided to seek approxim ; numerical solutions by 
finite differences. The governing equations are formulated 
in terms of the stream function, the vorticity and the po- 
tential function. The system of equations consists of two 
elliptic, three parabolic and an algebraic equation. The 
steady state solution is obtained by a large time solution 
of an unsteady problem. When a time dependent solution is 
sought, the solution procedure requires the solution of two 
elliptic equations after every time step. Since the solu- 
tion of elliptic equations is the most time consuming, it 
is important to select an efficient numerical scheme for 
elliptic equations. Both the successive-over-relaxation 
(SOR) and the alternating-direction- implicit (ADI) method 
are known to be efficient. However, ADI method of Peaceman 
and Rachford^ 2 ^ has been chosen because of its relatively 
high convergence rate compared to SOR as the number of nodes 

in the mesh becomes large. This has been shown in the 

( 8 ) 

numerical experiments of Birkhoff, Varga and Young. 

ADI is also used for the parabolic equations. The 
stationary ADI (fixed parameter) which is used throughout 
this study is proven unconditionally stable for a simple 
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diffusion equation by Peaceman and Rachford. Later this 

(27) 

property of unconditional stability is shown by Pearson 
to apply to the diffusion equation with first order deriva- 
tives included. Various numerical schemes used in various 
hydrodynamic problems are reviewed by Ames i ^ 

Once the steady state solution is obtained for the 
set of equations^ the pressure is calculated either by 
line integration of the momentum equation or by solving the 
Poisson equation obtained by taking the divergence of the 
momentum equation. The Poisson equation was solved by the 
method of Fourier analysis in which a computational 
algorithm known as fast Fourier transforms was employed. 

Summary of Differential Equations and Boundary Conditions 

The set of governing equations and boundary conditions 
given in Ch. II are summarized and rewritten below in a non- 
dimensional form? 

The stream function equation: 

) + i— (i— i|j ) s - i-rfH (47) 

l Pr r r 2'pr r z'z Z? 

a d 

The vorticity equation: 

H t - p [ H rr + f"r * * §!<«“«> r - §f <»''«> z] 

(continued) 


* The superscript (*) for the non-dimensional variables 
is suppressed from here on. 
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♦ 5£ U r » 

D 


[ e T- <u r + ?*k>] H + l!?[^! 


(r 2 ) 


Re 2 p, 


2 2 

1 Ro p 1 1 

(uu r + ^'^7 )+ ~T £(uw r + 5™z>J 


ap 


(49) 


The angular momentum equation: 


i r - ir + 

p rr r r ~ 2 1 zz 

t a 


!f<rpur) r 




Reer, 


(50) 


The mass diffusion equation: 


: t ° s c [ c rr + <? + iT> 0 r + pSz + p^z] 


Reer, 


•If? (rpuc) r - Bfe«> wc >z + ™‘ A 


bp 


(51) 


The potential function equation: 


♦rr + ( ? + r>*r + Vzz + 1 ?r*z + r 


P- r, 

4 

P 2 p 


0 


(48) 


,,^5? h'h ■'iii/\ii-iitii^lVn^ili ii infiii^Whiaii^, v-U *!i. .ill .w. . .'♦.. . »i... .. L. u ■. s *\ .', 
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The density equation: 


p AtP^ p B,P 

/Pn d + “ P A,P^ P B,P^ 


P A,P /P B,P 


The velocity equations: 


u “"IS^z " e *r 


w ■ Jpr^r * §♦* 


The pressure equations: 


P r * - F ! 


tp, - -P. 

az i 


V 2 P 


a -S 


where 


P 1 * p(uu r + ju 2 - 


Re t » 2 . 
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(52) 


(46a) 


(46b) 


(53a) 


(53b) 


(54a) 


h"> 2 » - 7> 


A. (u + li + ^ ) 

5Re ' r r srz'r 


(53c) 


p(uw r ♦ 2 w 2 > - |jV 2 w 


3aRe* u r + r + a'^z 


(53d) 


? (rP l>r 


+ ; (p 2>* 


(54b) 


And the boundary conditions: 


tj,(0,z) * 0 

(57) 

ifi(r,0) * 0 

(58) 

Z 

4,(l,z) * - / p(l,z) [-u in; j + <^ r (1/Z)] dz 

(59) 

4 i(r,l) = ij>(l,l) for b < r < 1 

(60a) 

♦ * “J / P Car r X) ar^X + |<|> 2 (r,l)] dr 


for 0 < r < b 

(60b) 

H (0,z) = 0 

X 

(61) 

H(r,0) " 'fc ( ‘ r** + *«’] ct.o) 


subject to * 2 - - [ji*Pr* r ] ( r, 0) 

(62) 
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Ha '” ’-l^rrj« 1( z, 



subject to 4» r | (1>2) • o 


(63) 

H(r ,i, . - |1 [<L* , ♦ 4-(- 

L K a pr 

P ^2 + ^22*] 

(r,l) 

subject to tj> 2 » - ^^|3pr$ r j 

(r,l) 

(64) 

r( 0 , 2 ) = 0 


(65) 

r(r,0) = 0 


(66) 

r(i,z) = i 


(67) 

r(r,l) a 0 


(68) 

c r ( 0 , 2 ) = 0 


(69) 

c 2 (r,0) = 0 


(70) 

ReS 

C r (l,2) + -£-- C -U in jC(l»2) a 0 


(71) 


c 2 (r,l) ■ 0 


(72) 
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♦ r (o,z) » a 

(73) 

<j> 2 (r,0) * 0 

(74) 

ij) ( 1 , z) a 0 

(75) 

<|> 2 (r,l) • 0 for b < r < 1 

(76a) 

<f>(r,l) * 0 for 0 < r < b 

(76b) 


Now that the parabolic equations are expressed in a general 
form the finite difference equations can also be expressed 
in a general manner. 

The parabolic equations in a conservation form can 
be written as: 

f t * Af rr + Bf r + Cf zz + Df z + E * r P uf ) r + F(pwf) s 

+ Gf + H (77) 

The coefficients A through H appearing in this equation can 
be obtained by comparison with individual equations. 

The Numerical Method by Finite Differences 

In order to apply the finite difference approxi- 
mations to the set of governing equations, the region of 
interest is divided into a finite number of subregions with 


O 


•V 




equal spacings. Since in the numerical computation the 

desired steady state solution is obtained by a large time 

solution of unsteady state problem , the time coordinate is 

also divided into a finite number of steps. With this grid 

system a differential equation of a continuous function is 

transformed into a system of finite difference equations of 

discrete mesh functions. The relationship between the 

continuous function f (r ,z, t) and the mesh function f? . is 

1 • j 

f i r j “ f < iAr » j**» nAt > 

since 

r^ « iAr, i * 0,1,2 ,,N 

Zj ° jA2, j B 0, 1, 2 ,.**... ,M 

t ft * nAt, n * 0,1,2, ,L 

Finite Difference Approximations of Parabolic 
Equations . The non-iterative implicit numerical scheme, 
the alternating-direction-implicit procedure of Peaceman and 
Rachford' is employed for the finite difference approx- 
imation. In this method, each equation is written in a 
sequence of two difference equations each of which is im- 
plicit in only one space coordinate involving only a half 
time step. Now the finite difference equations obtained 
after applying ADI to the general form of parabolic equations, 


Eq. 77 are; 

for the first time step. 
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2,-n+l/2 -n . A,-n+l/2 --n+1/2 ^+1/2. 

a' r i,j E i,j» ’ 2f i,j *1+1, l’ 


4. ^ • gn/fH Jf*' + f® } 

+ 2 lf i+l,j *i-l, j' + PCl *i # j-l 2f i,j + f i, j+1 


hD,-n 


(f n - f 11 > + £If r 0 n+l/2 n+l/2-n+l/2 

U i,j +1 r i,j-r + 2 lr i+l p i+l,j u i+l,j f i+l,j 


' r i-l°i-l(j U i-l{j f i-l{j ) + F“ <P i, j+l w i, j+l f i, j+1 

+ 3 a<f itj /2 + *;.!> + « < 78a > 


and for the second half step. 


2, -n+1 -n+1/2. A.-n+l/2 --n+1/2 . -n+1/2 

■ f i,j >■ * 2 f i,j + *i+i,j 

♦ r lf HM - + * c(f !ft-i - 2 f ?!i + f iti+i> 

♦ ?ku - *&. i> + ^ i+i«<ft*i:ft 

„ n+1/2 n+1/2 -n+1/2. . hP, n+1 n+1 -n+1 

" r i-l p i-l,j U i-l,j f i-l,j' + 2~ (p i,j+l W i,j+l f i,j+l 

n+l n+1 -n+1 . . 6G,-n+l . -n+1/2. . » u naM 

♦ Pi ( .i-i w i,j-i f i,j-i ) + + £ i,j > + SH (78b) 




nTY p™.-: T l. n,,;^^,-, tr ,- TV - ,5,-,-f,^: TT , TI - YtV „- rT , J 



where h * Ar, k * Az, a = At/hk, @ = h/k and 6 » hk. 

Above equations are rearranged and given in a form conve- 
nient for numerical computation: 


+ h U f T,Y 2 * °i,i *Ki3 • d i.j < 79a) 


a** f n+ i + k** f n +l . -** p n+l 
a i,j f X,j-X + + i, j £ i, j+X 


d i:i 


(79b) 


_* = la kjj k_ n+1/2 n+1/2 

a !,j 5* 5® 5® r i-x p i-x,j u i-x,j 


b I.J 


.2 2, . o. 

o 6* + 5° 


c J,j " I* + 5 B + 5® 


(80a) 


d Lx 


- -f? . - (f? . . - 2f? . + f” ...jec 

a i/j i»j-l i#D i*j+l 


«h/*n *n , n h , n ,,n -n 

7 (f i,j+l f i» j-1 D " ? ip i,j+l w i,j+l f i, j+1 


6^, -n 


- **i.j ‘ 6H 


£ P n+1 w n+l 

I Pp i,j-l w i,j-l 


(continued) 






b Hj 


- | - 20C + ~G 


C ** 

°i,j 


ec + So + fe , P? + L,w? + f J 


2 if j+1 i# j+1 


(80b) 


<j** - - 2^f*' + l/2 _ l,-n+l/2 2 f n+ ^/2 + f n+ ^/^)A 

a ifj " a r ifj $ U i-lfj ^ £ ifj + *i+l,j ,A 


. k (f .n+l/2 _ -n+l/2. fi _ k, n+ 1/2 n+l/2-n+ 1/2 
2 u i+lfj f i-lfj ,B T tr i+l p i+l,j tt i+l, j £ i+lf j 


- - Kj /2 - SH 


Finite Difference Approximations of Elliptic Equations* 


The ADI of Peaceman and Rachford is also used for the ellip- 


tic equations. Iterations are carried out first implicitly 


in r coordinate and next in z coordinate in an alternating 


manner. Following Birkhoff , Varga and Young 1 ' for the case 


of a single iteration parameter , the finite difference 


equations can be expressed as the following. 


The finite difference equations of the stream function: 


for the first iterate , 


IP* 4*4 - P 


r r *?_, 4 ♦ lh~ z- 1 


i*j i i 




j r i+l + p i,j r i 


^fj^ + °i-lf j r i-l >+ q ®] * i,j 


[p i+l,j r i+l + p i # j r i > 




- i i ii.' MM P' W * I ' l ii h i. » jf w— wpmpwy 
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♦i+l,j ° q s*i,j + 


16 r 1 

a r A l p i,j+l + p i,: 


( *i,j+l * W 


p i 


' *1.1-1*] + ^*1*1.1 


(81a) 


and for the second iterate , 


2S __ 

tt* 

*i.j-i 

a2,r i (P i,j + ^j-l* 

1 1 

** 

•i.j + •i.i-i > * 9 *J 

*i,i ■ 

_ * . 2 r 



+ e [ p i+l, j r i+l + p i , j r i 

.?l +V-T7I rS'i * + 


(81b) 


The finite difference equations of the potential function: 
for the first iterate. 


[- > • M; • -‘“ii.;,; 1 - 1 ' 1 '] • <t • y.;,. 

■ [s * I l r7 * *(♦!,) ' Vl.) ‘ 


f 


2 *i,j + + 

XiJ 

p i,j 


p i,j-l> ( *i,j+l “ <|, i,j-l > 


(82a) 


and for the second iterate , 


p i,i+l " Pijj-l 


+ v*Hj ' t ?* 1 


p iti+i "filial 


* vi.j + iN-M ■ *i.i 


+ +i+i, j 1 + !% ♦ <+i+i,j - 

+ 

p i,j 


fiiiw,.* 


(82b) 


In above equations the functions with superscripts (*) and 
(**) represent the first and the second iterate respectively. 


The Numerical Method and the Finite Difference 


:oxi- 


m at ions of the Boundary Conditions , The finite difference 
equations obtained by ADI method as given above result in a 
tridiagonal system of algebraic equations since each equation 
is implicit in only one direction. The solution of the tri- 
diagonal system is obtained by an efficient matrix inversion 


v . - • 


known as the Thomas method described by LeeP' L; In general, 
non-linear problems can be treated as a set of linear prob- 
lems. In the present problem non-linearities exist in the 
equations as well as in the boundary conditions, i.e., the 
boundary conditions for the vorticity. When it is necessary 
to account for non-linearities, the equations and the bound- 
ary conditions are iterated over the same time step until a 
predetermined convergence criterion is satisfied before ad- 
vancing to a new time level. 

The boundary conditions are all approximated by finite 
differences with second order accuracy to be consistent with 
the accuracy associated with the finite difference equations. 
The first and second order derivatives for boundary points 
are expressed in finite differences as in the following. 

The first order derivative* 


f x (x o } = [** 3f(x o ) + 4f(x o + Ax) “ f(x o + 2Ax) ] 

+ 0 (Ax 2 ) 

f x (x o ) = lb? [ 3f(x o ) ~ 4f(x o " Ax) + f(x o * 2Ax) ] 


+ 0 (Ax 2 ) 


The second order derivative with a specified value of first 
order derivative: 




£ xx <x o > ' ~2 [-6f(x Q )Ax - 7f(x 0 ) + 8f(x 0 + Ax) 

-f<x o +2Ax)] + 0(Ax 2 ) 

f xx* x o) " £x' ( x o )ax - 7f (x Q ) + 8f (x Q - Ax) 

2^x 

- f(x o - 2ax) ] + 0 (Ax 2 ) 

Finally, the second order derivative with no restriction 
on the £irst order derivative: 

f xx <x+> “ "T* [ 2f ( x Q ) ** 5f < x o + Ax) + 4f (x Q + 2 Ax) 

- f (x Q + 3ax) ] + 0 (Ax 2 ) 

f xx < x ;> ° "T* [ 2 f(x Q ) “ 5f(x o - Ax) + 4f(x - 2 Ax) 
Ax 

- f(x Q - 3AX)] + 0 (Ax 2 ) 


In the following, the computational procedure is de- 
scribed to obtain the near-steady-state solution by the time 
dependent formulation. At a certain time level, n, all the 
functional values are known from which functional values at 
a new time level, n+1, are to be obtained. This procedure 
is repeated until a predetermined criterion is met for the 
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steady state solution. 

(1) Compute r solving the ADI equation of r using 
the functional values known at time level n. 

(2) Compute H solving the ADI equation of H with 
boundary conditions based on current value of 

ip and p . 

(3) Using the elliptic ADI equation of ip compute 

by repeating the iterations* until the following 
convergence criterion is satisfied: 



where e.^ is some small positive number and ip . i 

3TG f 

is some reference value of +h e s tre am funct ion. 

(4) Improve u n+1 and w n+1 using current value of ip j 

computed in step (3) . j 

(-S4— Gempute c n+1 solving the ADI equation of c using j 

the currently available values of other functions 

involved. 

n+l 

(6) Update p using recent value of c obtained in j 

step (5). j 

(7) Using the elliptic ADI equation of <p, compute $ j 

by repeating the iterations* until the following j 

| 

* If only the steady state solution is needed, the iteration 
is terminated after one cycle, which yields a meaningless in- 
termediate solution. However, it results in an appreciable 
reduction in computation time. 

I 

1 


9 
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criterion is met: 

|*"!j * ♦!,j| 5 e 2 *ref 

where is some small positive number and 6 - is 

* ref 

some reference value of the potential function. 

( 8 ) Update u n+ ^ and w n+ ^ using current values of \p , p 
and $ obtained in steps (3) , ( 6 ) and (7) respec- 
tively. 

(9) Evaluate the boundary values of H based on current 
values of p and 4 . 

(10) Return to step ( 1 ) and repeat the cycle without 
advancing the time step*. With these non-linear 
time iterations the solution approaches that of 
the time dependent problem. The iteration is term- 
inated when the following convergence ' criterion is 
satisfied. 

u av I 1 k + 1 .k I < 

Max. |» 1(j - * 1#j | - c 3'l , r ef 

where k is the iteration counter and e 3 is some 
small positive number. This completes the com- 
putational procedure for one time step advancing. 

(11) For a steady state solution, above procedures are 
repeated until the time dependent solution convezges 

* If interested in large time solution only, this iteration 
is not necessary except in the first few time steps. 


r 


to the steady state solution. The convergence 
criterion is again * 


Max. 




< 

■ e 


4* ref 


where is some small positive number. 

Pressure Calculation 

Once the steady state solution is obtained from the 
set of governing equations formulated in the stream function , 
vorticity and potential function , the pressure is obtained 
either by a line integration or solving the Poisson equation 
of pressure. 

Pressure by Line Integration . This method involves 
evaluation of pressure gradients at every grid point in the 
region and numerical integration along conveniently chosen 
paths. In this study, pressure is firs±_in£egrated axially 
along r* = 0.5 from z* * 0 to z* « l f and next radially 
inward and outward to cover the entire region. Pressure 
gradients are obtained using Eq. 53 for known velocity 
field. This equation is rewritten below: 

He 2 2 

P r - - P(uu r ♦ jk<V*u - Hj) + ^ 


(u + — + ) 

' r r aVr 


(53a) 
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> z - -P<uw r + + |jV 2 w 


1355 <u r 


+ H 
r 


+ ~w 


Z^ z 


(53b) 


The integration formula used here is known as Euler's 
polygonal curves: 


f (x Q - Ax) - f(x 0 ) + £ [-Sf x (x o - Ax) - 8f x (x 0 ) 
+ f x (x G + Ax)] 

f(x 0 ) . f(x 0 - Ax) + ff [*£ x (x 0 - Ax) + 8f x (x 0 ) 
-S x (x 0 + ax)] 


£(x 0 + Ax) * f(x G ) + jj [" f x < x o * Ax) + ® £ x (x o ) 


+ 5f x (x o * 4x) ] 

Pressure by Solution of Poisson Equation . The pressure 
can also be obtained by solving the Poisson equation of 
pressure Eq. 54. This equation is rewritten below: 

„2 

V p » -S(r,z) ( 54a ) 


where the source function 8 is a known function when the 


/ 

l 


t 
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velocity field is obtained. The boundary conditions are 
obtained from the known velocity field by Eq. 53a and 53b. 

Since it is difficult to solve a Poisson equation with 
Neumann boundary conditions all around , it was decided to 
obtain solutions by a series expansion in eigenfunctions. 

This method is similar to that of Green's function and is 
chosen because of its computational efficiency. 

The Poisson equation with inhomogeneov. s boundary con- 
ditions is first converted into one with homogeneous boundary 
conditions. This is accomplished by constructing a simple 
function p st’/ih that it satisfies the original inhomogeneous 
boundary conditions f 18) Then the solution to the homogeneous 
problem p is related with that of the original inhomogeneous 
problem p by: 

P = P ” P D (84) 

The resulting Poisson equation with homogeneous Neumann 
boundary conditions is: 

V 2 p » -S (85a) 

where 

S = S + V 2 p^ (85b) 

* o 

Now assume a form of the solution as a finite sum of Fourier 
cosine functions of z and at the same time expand the source 
function in a finite cosine series: 


S 


I 


( 86 ) 


P 


K <r > + 


N 

l 

n=l 


c (r)cosntfz 
n 


— 1 w 

“S a ^d (r) + l d (r)cosnirz (87) 

n=l n 

The number of terms in the series is taken same as that of 
the mesh points in z coordinates.^ 5 ^ Substituting these 

equations into Eq. 85a yields a set of ordinary differential 
equations for the Fourier amplitudes : 

i 2 2 

c n (r) + r°n (r) “ ^Hr c n (r) 88 d n (r) ' n “ 0,1,...., N 

a 

( 88 ) 


with homogeneous Neumann boundary conditions: 

c^(r) « 0 at r * 0 and r = 1 (89) 

This set of ordinary differential equations was solved 
numerically by finite differences, in which the Thomas method 
is used for the resulting tri-diagonal system. A matter 
deserving a comment here is that the fundamental Fourier 
amplitude c Q (r) can best be determined within an additive 
constant. This is the reason why the pressure can only be 
determined within some constant. Once p is obtained, p is 
determined by Eq. 84. 
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For purposes of computational efficiency, a simple ver- 
sion of the fast Fourier transforms (FFT) is used. The 
reduction in computation time is accomplished by using the 
symmetry property of cosine functions to minimize the multi- 
plicative operations. The reduction in computation time 
realized is approximately by a factor of 4. The method of 
performing the fast Fourier cosine transforms is given in 
Appendix A. 


CHAPTER IV 


ERROR ANALYSIS 

In this chapter, various errors involved in the 
numerical computation are discussed. Efforts are made to 
check the solutions obtained internally as well as extern- 
ally by way of comparisons with other related studies. Un- 
fortunately, a solution to the present problem in its full 
complexity does not exist in the published literature. How- 
ever, results of visualization experiments and numerical 
studies of similar rotating flows are found and they are 
used for both qualitative and quantitative verifications. 

There are two major sources of error in solutions 
obtained by numerical methods. The discretization error 
is the difference between the exact solutions of the finite 
difference equations and those of the differential equations. 
In addition, actual numerical solutions include another 
error called the round-off error since the computing machine 
uses only a finite number of digits. 

In this work, equations are posed as an initial value 
problem and steady state solutions are obtained by contin- 
uing the computation until the solution no longer changes. 

The solutions thus obtained contain both discretization and 
round-off errors. If, in the solution procedures, these 
errors remain bounded as time increases, the numerical 
scheme is said to be stable. As the mesh size decreases, 
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the finite difference equations approach the differential 

equations since truncated Taylor series expansions are used 

for the finite difference representation. The consistency 
( 29 ) 

condition is thus satisfied. For certain simple linear 
problems , the consistency and stability are shown by Lax's 
Theorem' to imply convergence. Hence for simple prob- 
lems, it is proved that stable numerical solutions converge 
to the solutions of differential equations when the mesh 
size approaches zero. However, the theorem does not extend 
in general to apply to more complex and non-linear equations 
Therefore, the result of Lax's Theorem can only be used as 
guides for non-linear problems and the numerical solutions 
obtained in this work are subject to verification as to 
their validity as solutions to the differential equations. 

In order to answer the question of how closely the 
numerical results obtained approximate the true solutions 
of the differential equations, a series of checks are made. 
They a e briefly listed below and described in the remainder 
of this chapter. 

(1) Verify how well the numerical solution satisfies 
the constraint of the continuity equation on a 
discrete basis. 

(2) Check how closely the compatibility condition* 
is satisfied by individual species (the over-all 

* The surface integral of flux on boundary being compatible 
with the volume integral of source in the interior. 
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species balance) . 

(3) Check how well the numerical solution satisfies 
the compatibility condition for the Poisson 
equation of pressure. 

(4) Establish the uniqueness of the steady state 
solution by using different initial conditions. 

(5) Compare the numerical results of the present 
study with reported visualization experiments 
for a qualitative verification. 

(6) Use the computer program of this study to solve 

(25) 

the rotating flow problem of Pao and compare 
the results for a quantitative verification. 

(7) Estimate the discretization error by an extra- 
polation using different grid sizes. 

(8) Calculate the pressure by both line integration 
of the momentum equation and by the solution of 
the Poisson equation for the pressure and com- 
pare the results. 

The first three of these checks are considered to 
be weak internal checks. Whereas the next four are strong 
positive checks. Finally, item (8) is believed suffrcient- 
for the verification of the pressure results. 

Internal Checks 

The checks mentioned in this section are simple in- 
ternal checks which partially contribute to the verification 
of the numerical solutions reported in this work. The first 
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verification is of the continuity equation in the finite 
difference approximation. The error in the continuity 
equation is: 

e (r,z) * I(r p u) r + |w z - er A (90) 

where e (r, z) represents the error in a non-dimensional form. 
The relative local error, i.e. the absolute error e(r,z) 
divided by the largest of the three terms involved, are 
computed at every interior grid point for four selected 
runs . These four runs were selected in such a way that 
they are representative of all the cases considered in this 
study and also they are expected to have the maximum error. 
The results show that the relative errors are less than 
10 -5 at all interior node points except at the two nodes 
adjacent to the boundary point at which the exit hole ter- 
minates. At these two interior nodes the errors are found 
to be as large as 50 percent. It indicates that the viol- 
ation of the continuity constraint is localized very close 
to a point which is nearly a singular point and does not 
extend into the interior. 

The second check made is on the over-all balance 
of species. For every computer run, the total amount of 
interior mass injection and the sum of mass transfer across 
the boundary are computed for the component A. This checks 
the compatibility condition for the mass diffusion equation 
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as well as the general validity of the numerical solution. 

results of this check — shows that the condition is sat— 
isfied to within less than 5 percent error. 

The third simple check carried out is on the compati- 
bility condition of the Poisson equation of pressure. The 
source term and the boundary conditions are evaluated from 
known solutions of kinematics. If the numerical solution 
obtained for the kinematics is valid , the bo un dary conditions 
are expected to be compatible with the differential equation. 
This is the reason why the verification of the compatibility 
condition serves as a check for a valid solution of kine- 
matics. The result shows that this condition is met within 
less than 10 percent deviation for all rotating flows con- 
sidered. For non-rotating flows, however, the error is found 
to be as large as 50 percent. For this reason, the line 
integration method was adopted for pressure calculations 
of all non- rotating flows. 

Finally, the uniqueness of the steady state solution 
is established since identical results are obtained from 
three different initial conditions for two chosen flow 
cases. The initial conditions considered are those of 
solid body rotation, those in which the interior values are 
interpolated from boundary conditions, and finally initial 
values of zero everywhere in the region. Single component 
flows of Re * 20, Re t « 0 and Re * 20, Re^ ■ 50 were used 
for this part of the analysis. 
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Comparison with Results of Prior Works 

In this section comparisons are made of the results 

of the present work with those of prior studies of other 

investigators. The first comparison is made with the ex- 

(30) (32) 

perimental visualization results by Roschke and Ross. 
These studies were concerned with the secondary flow pat- 
terms obtained in jet-driven vortices which had a small 
exit hole on one end. Schematic sketches of the over-all 
flow pattern were drawn from the visualization results. 

These sketches showed features such as boundary layer flows 
on stationary end walls and vortex core flows along the 
cylinder axis toward the exhaust hole. This over-all second- 
ary flow pattern is consistent with the numerical results 
of the present study as given in the next chapter * And 
this is considered to be a global, qualitative verification. 

The second check made is with the numerical solution 
of Pao^ 25 ^ who considered the steady motion of an incom- 
pressible flow in a finite cylindrical container. In this 
problem one end wall was stationary and the other rotated 
with the periphery. The governing equations of this prob- 
lem are identical to those of the present study for the 
special case of single component flow. Furthermore, this 
problem is considered to possess the same mathematical com- 
plexity as well as the flow features of the present problem, 
the boundary conditions of Pao's problem are: 
r = 0; u = 0, v = 0, w f = 0 




r » l; u = 0, v = 1, w = 0 


z=0;u=0, v=0, w=0 
z * 1; u«0, v»r, w*0 

The computer program written for the present problem was 
modified to meet the above boundary conditions. This mod- 
ified program was used to solve Pao's problem for the case 
of Re fc = 100* and the result obtained was compared with that 
given by Pao. In Figure 2 the stream function is compared 
for a fixed axial location , z* * 0.5. The axial velocity 
along the cylinder axis is compared in Fig. 3. Finally, 
a comparison is made in Fig. 4 of the radial velocity a- 
long three different radial distances, r*®0.2,0.4 and0.6. 

The comparison is excellent as shown. It may be 
mentioned that in both Pao's and the present work, the fi- 
nite difference approximations were of second order accuracy 
although different numerical schemes were used. Pao's 
numerical results are supported by his visualization ex- 
periments. The comparison with Pao’s solution is considered 
to serve as a strong quantitative verification on the valid- 
ity of the numerical results reported in this work. 

Estimate of Discretization Error 


In this section, discretization error is estimated 
based on a quadratic extrapolation which is similar to 

/ 28) 

Richardson's extrapolation. The method involves 


* The characteristic velocity and the length used in the 
tangential Reynolds number are v c and R respectively. 
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calculating the approximate mesh function f(h) for three 
different mesh sizes. From these calculations the exact 
solution f(0) corresponding to zero mesh size and also the 
discretization error Eg can be evaluated. Suppose the 
mesh function f which depends on grid size is calculated 
for three different grid sizes, h, eh and 6h where 
0 <6 < e <1. Then from Taylor series expansion, 

f (h) « f(0) + ah + bh 2 + 0(h 3 ) 

f (eh) = f (0) + aeh + be 2 h 2 + 0(h 3 ) (91) 

f (<$h) * f (0) + a<$h + b6 2 h 2 + 0(h 3 ) 


where a and b are constants which depend on the first and 
second order derivations of f respectively. But in this 
analysis the order of approximation in all the finite 
difference equations is of 0(h) and hence a * 0. There- 
fore the ratio, 


f (h) - f (eh) m (1 - e 2 )bh 2 ± 0(h 3 ) 
f(eh) - f(6h) ( e 2 - <s 2 )bh 2 + 0(h 3 ) 


(92) 


approaches the limit as the mesh size h tends to zero: 


f (h) - f (eh) 1 - e 2 

f (eh) - f(6h) * “5 

C — 0 


( 93 ) 


The error Eg is obtained from Eq. 91 and 92: 
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E d « f(h> - f(0) « f (94) 

G “6 

This method of error analysis is applied to a non-rotating 
flow of Re = 20 using four different grid sizes, 1/20, 1/28, 
1/36, 1/44. From the results then the errors are estimated 
for grid sizes of 1/20 and 1/28. The errors were calculated 
based on the computed values of circulation along six 
chosen paths as in the following t 

1 

/ w*dz* along r* * 0.25, 0.50 and 0.75 

0 

1 

-/ u*dr* along z* = 0.25, 0.50 and 0.75 

o 

The integral quantities were chosen rather than local 
quantities since they reflect errors associated with a 
group of mesh points along the path of integration. 

The results show that the relative discretization 
errors associated with grid sizes of 1/20 and 1/28 are 12 
and 4% respectively. For rotating flows, however, the 
errors are expected to be higher because of the induced 
secondary flows and thin boundary layers which have steep 
gradients of velocity. Nevertheless, for the range of 
flow parameters considered in this study it is believed 
that the errors are not much affected by the wall rotation. 
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The calculated circulations along the chosen paths and their 
variation with mesh spacing are shown In Pig. 5 and 6. 

Check on Pressure Calculation 

Pressure is calculated both by line integration 
using the integration formula, Eq. 83, and by solving the 
Poisson equation of pressure. A comparison of the results 
is given in Fig. 7 for the case of Re * 1 and Re^ = 50 of 
the single component flow. In this figure the pressure is 
plotted against radial distance for three chosen axial 
locations. The results show good agreement between the 
two methods and this is considered sufficient to establish 
the validity of the results for the pressure calculations. 
Concluding Remarks 

In light of the results of the error analysis that 
are discussed in this chapter, the numerical results re- 
ported in the next chapter can be viewed with confidence. 

As pointed out in the earlier section of this chapter, the 

details of flow solutions in the vicinity of the exhaust 
hole are not considered to be accurately discerned. How- 
ever, the associated inaccuracy is very much localized and 
limited to the immediate vicinity of the exit port. This 
was brought out by the result of the check on continuity 
in the previous section. Finally, it is shown by the re- 
sults of the compatibility check for the Poisson equation 
of pressure that the line integration method should be used 
for all non-rotating flows, whereas for rotating flows, that 
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CHAPTER V 


NUMERICAL RESULTS 

Results of the numerical calculations are presented 
in this chanter. The computer program that was used to 
produce them was written in Fortran V for the UNIVAC 1108 
and included in Appendix B. In early stages of data ac- 
quisition, the limits of stable solutions were explored for 
the governing parameters, particularly the axial and tan- 
gential Reynolds numbers. The result showed that stable 
solutions were obtainable without excessive machine time 
for Reynolds numbers upto 20 for rotating flows with tangen- 
tial Reynolds numbers as large as 100. Whereas for non-ro- 
tating flows stable solutions were obtainable for Reynolds 
numbers as niqh as 2000. 

The maximum time step and the optimum iteration 
parameters for given flow cases are determined by numerical 
experiments for the parabolic and elliptic equations re- 
spectively. The exhaust hole radius is kept constant 
r e * 0.125R for all the flow cases considered. Aspect 
ratios considered are 1 and 2 and the Schmidt number is 
fixed, i.e. S c * 1 for all binary flows studied. 

Although it is possible in principle to compute the 
required mass addition rate of component A that is necessary 
to maintain a predetermined critical mass in the chamber, 
the computation procedure required to do this is difficult. 
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Therefore, for purposes of simplicity the containment of 
component A is calculated for given values of the mass 
addition rate. Since it is desirable to keep the mass 
addition of simulated fissionable material as small as 
possible, relatively small values, 0.01 and 0.05, are con- 
sidered for the injection rates of component A(f A _ B ) . 

Grid systems of 29 x 29 and 21 x 41 are used for 
flows of aspect ratio of 1 and 2 respectively. In Table 1, 
single component flows that were considered in this study 
are summarized. Whereas Table 2 lists the binary fluid 
flow cases studied. The component simulating fissionable 
material is introduced by a uniform interior source dis- 
tributed over an annular containment region, 0.125 < r* < 
0.875 and 0.125 z* < 0.875 for all binary flow cases. 

This choice of mass injection schedule is made in order to 
reflect local effects of secondary flows on the resulting 
containment factor. Although the assumed uniform source 
distribution may not be an efficient one for achieving high 
containment, it is considered most suitable for studying 
the effects of local secondary flows on the retention of 
fissionable material. 

Because only the steady state solution is of interest 
it was necessary to iterate the elliptic equations only 
once for each time step. Pinal steady state solutions are 
obtained with a minimum of computer time in this way, but 
the results produced before convergence is reached does 
not represent the true transient behavior of the flow. 



Table 2 . Summary of Numerical Results for Binary Fluid Flows . 





Results are presented in the following three sections 
separately for single component flows , binary fluid flows 
and for the pressure respectively* Calculated values are 
dotted along selected r* and/or 2 * values whichever is 
deemed appropriate for the particular case under consider- 
ation. Results are given in a manner that they show para- 
metrically the effects of rotation on the resulting second- 
ary flow, concentration field and finally the containment 
factor. 

Single Component Flows 

Plows with an aspect ratio of 1 are taken up first 
for discussion. Streamlines for flows of Re = 1 and for 
Re t • 0 and 50, are presented in Fig. 8 and 9 in which the 
overall flow configurations can be seen. With no rotation, 
i.e. Re t = 0, the fluid enters at the porous wall at the 
periphery and flows more or less uniformly throuah the bulk 
of the chamber. However, with rotation, Re^ * 50 , a pair 
of recirculating closed flow cells are formed in the main 
body of the chamber. The primary fluid injected at the 
periphery flows in the boundary layers on the stationary end 
walls and around the recirculatina closed cells. The mass 
flow rate in the secondary recirculating cells is measured 
by the quantity, |ij/| ma x the closed cells. This quantity 
increases from 0 to 45.9 as the wall rotation, Re^ varies 
from 0 to 100 for fixed primary flow rate of Re ■ L This 
is shown in Fig. 13. 

It is to be noted that for a particular ratio of 
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wall rotation to primary flow, i.e. Re t /R e ■ 50, a closed 
streamline appears, which contains the chamber center 
line. The appearance of this particular streamline brings 
about a new flow feature worthy of discussion* Half of 
the primary flow that flows through the left end wall 
boundary layer, changes direction, passes through the 
central bulk region and joins with the remaining half to 
finally exit through the exhaust hole. This flow con- 
figuration is not a desirable one for fluid dynamic con- 
tainment. It may also be mentioned that this flow can not 
be studied by the simple analytical model of Rosenzweig, 
Lewellen and Rosst 31 ) since they assumed the axial veloc- 
ity to be uniform in the region of the central vortex core. 

Streamline patterns for a higher Drimary flow rate, 

R e * 20, for wall rotations of Ret«0, 50 and 100 are 
given in Fig. 10, 11 and 12 respectively. In this case a 
similar partern of vortex motion is set up by the wall 
rotation. The primary flow is through the end wall bound- 
ary layers and around the outside of the induced recir- 
culating flow cells. The secondary flow rate, |<H max of the 
closed cells varied from 0 to 2.64 as the tangential 
Reynolds number increased from 0 to 100. This is shown in 
Fig. 13. A comparison of this result with that of R e * 1 
indicates that the secondary flow rate decreases for a fixed 
rate of wall rotation as the axial Reynolds number increases. 
The streamline patterns show that the fraction of primary 
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flow that passes through the central bulk region decreases 
when the tangential Reynolds number increases from 50 to 
100. It may be stated from the streamline patterns for 
both R e * 1 and 20 that increasing the tangential Reynolds 
number will at first decrease the fractional primary flow 
turning to the main body when the axial Reynolds number is 
fixed. This tendency persists until a minimum is reached. 
Further increase in rotation will increase the fractional 
primary flow until a flow condition develops where the 
closed streamline appears as given in Fig. 9. This sug- 
gests the existence of an optimum rotation rate for a 
given primary flow that yields the most favorable flow 
structure for containment. 

In the case of an axial Reynolds number of 20 , a 
small closed flow cell that contains the chamber 
axis is formed in front of the exit port. This has not 
been observed in the case of Re * 1. A close examination 
of the pressure distribution along the axis revealed a 
high pressure point immediately downstream from the cell. 
It is thought that the porous plua conditions imposed on 
the orimary flow at the exhaust hole causes a region of 
local high pressure on the center line so that the primary 
flow can adjust within a short distance to comply with the 
boundary conditions. The local high pressure induced by 
the exit port boundary xmditions then causes a local re- 
verse flow to form a closed flow cell. In the case of 


Re • 1, however, the pressure level required for the neces- 
sary adjustment was not high enough to cause a back flow. 
However , the detailed flow configuration as given in the 
vicinity of the exahust hole should be viewed with reser- 
vations. The grid system used in the numerical computation 
yields only three and a half nodes across the exit port 
and tli us the rapid changes in the primary flow in this lo- 
cal region have to be approximated by only a few points. A 
local high value of the error in the continuity constraint 
was also observed in this region as discussed in Ch. IV. 

Radial and axial velocity components are plotted for 
Re * 1 in Fig. 14 and 15 for Re t * 0 and 50 respectively. 
These figures support the discussions given above in regard 
to the streamline patterns, with no rotation, an asymmetry 
in the veloicty distributions about the mid-olane, z* * o,5 
is shown. At Re t « 50, the relative effect of the primary 
flow is low and as a result velocity distributions are near- 
ly symmetric about the mid-plane. The radial velocity dis- 
tribution shown in Pig. 1.5 reveals that the end wall bound- 
ary layer thickness grows rapidly as the flow approaches 
the chamber axis. The results of Re t • 100 with Re » 1 
are not given here since they are similar to the case of 
^®t a but with an increased intensity of the recirculating 
flow features described above. 

The distribution of angular momentum is shown for 
Re = 1 and Re t * 50 in Fig. 16. It indicates a radial 


distribution which does not vary axially in the main oody 
of. the chamber . Hence a steep axial gradient exists near 
the end walls. This supports the T(r) assumption made by 
Rosenzweig, Lewellen and Ross< 31) in their analytical model 
of driven vortex flows. The observed high axial gradient 
of angular momentum near the end walls provides a large so- 
urce term in the vorticity equation. Through the stream fu- 
nction equation, this leads in turn to a large variation in 
the vorticity and results in boundary layers in this region. 
The distribution of angular momentum also is shown to be 
symmetric about the mid-plane r . 

Plots of radial and axial velocities for Re * 20 
are given for Re^. * 0, 50 and 100 in Pig. 17, 18 and 19. 

The velocity distributions are similar to those for the case 
of Re ■ 1. But the distribution is less symmetric about 
the mid-plane because of the higher primary flow rate. 

This effect is more pronounced for lower tangential Reynolds 
numbers and in the region close to the axis. The angular 
momentum is given in Pig. 20 and 21 for Re t * 50 and 100 
respectively. These figures show distributions similar to 
the case of Re * 1. The symmetric distribution of angular 
momentum about the mid-plane is not altered by the increased 
orimary flow rate. 

In the rest of this section results for an aspect 
ratio of 2 are discussed. Only a few runs are made for 
this case . The streamline patterns of Re*l are shown in 
Fig. 22 and 23 for tangential Reynolds 
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numbers of 0 and 15 respectively. For the particular value 
of the ratio, Re fc /Re ■ 15, the pattern of a closed stream- 
line that joins with the chamber axis is aqain formed and 
it causes a larqe fraction of the primary flow to pass 
through the central reqion of the chamber. The streamlines 
for Re « 20 are given in Fig. 24 and 25 for cases of Re t = 0 
and 15 respectively. For these values of parameters, the 
fractional primary flow turning to the central region is re- 
duced because of the increased axial Reynolds number for 
the fixed tengential Reynolds number of 15. The secondary 
flow rate for this case in which ;he aspect ratio is 2 is 
also plotted in Fig. 13, From this figure, it may be seen 
that the secondary flow rate increases as the aspect ratio 
increases for fixed axial and tengential Reynolds numbers . 

Velocities are plotted for Re * 1 in Fig. 26 and 27 
for Re^. - 0 and 15. The velocities for Re * 20 are given 
in Fig. 28 and 29 respectively for Re t « 0 and 15. The 
velocity distribution is similar to the case with an aspect 
ratio of 1. The asymmetrical distribution of radial and 
axial velocities can be seen for the case of Re = 20 and 
Re fc =15. It appears that thinner boundary layers are 
formed on the end walls. Angular momentum is given in Fig. 
30 and 31 for Re = 1 and 20 respectively. The reqion in 
which the axial gradient of angular ...omentum is steep also 
seems narrower than that for the case of aspect ratio of 1. 
This results in a higher axial gradient of angular momentum 
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and hence a higher secondary flow rate. This is in agree- 
ment with the result given in Fig. 13. Finally, the general 
distribution of angular momentum is similar to the case 
with an aspect ratio of 1. 


Axial distance, z 


Fig. 9. Streamline Pattern for Re = 1, Re.= 50 























•Mill 


Axial distance, z 


£L%^&2sttCes'±gC}G**‘ -' 3=0 


Radial distance# r 

Fig. 14, (a) Radial Velocity vs. Axial Distance for 
Three Selected Radial Locations, (b) Axial 
Velocity vs. Radial Distance for Three Select 
ed Axial Locations, for Re — 1, ® and 
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Pig. 15. (a) Radial Velocity vs. Axial Distance for 
Three Selected Radial Locations, (b) Axial 
Velocity vs. Radial Distance for Three Sele 
cted Axial Locations, for Re 8 1, Re fc * 50 
and a 8 1. 
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Fig. 17. (a) Radial Velocity vs. Axial Distance for 
Three Selected Radial Locations, (b) Axial 
Velocity vs. Radial Distance for Three Sele 
cted Axial Locations, for Re * 20, Re t * 0 
and a * 1. 
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Fig. 18. (a) Radial Velocity vs. Axial Distance for 
Three Selected Radial Locations, (b) Axial 
Velocity vs. Radial Distance for Three Sele 


cted Axial Locations, for Re 
and a * 1. 
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Pig. 19. (a) Radial Velocity vs. Axial Distance for Three 
Selected Radial Locations# (b) Axial Velocity 
vs. Radial Distance for Three Selected Axial Lo- 
cations. for Re * 20# Re* * 100 and a • 1. 
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Fig. 26. (a) Radial Velocity vs. Axial Distance for 
Three Selected Radial Locations, (b) Axial 
Velocity vs. Radial Distance for Three Selected 
Axial Locations, for Re « 1, Re. ** 0 and a = 2. 
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Fig. 27. (a) Radial Velocity vs. Axial Distance for 
Three Selected Radial Locations, (b) Axial 
Velocity vs. Radial Distance for Three Selected 
Axial Locations, for Re * 1, Re. = 15 and a « 2 
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Fig. 29. (a) Radial Velocity vs. Axial Distance for Three 
Selected Radial Locations, (b) Axial Velocity 
vs. Radial Distance for Three Selected Axial Lo- 
cations, for Re = 20, Re. = 15 and a = 2. 
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Binary Fluid Flows 

Results for binary fluid flows are presented and 
discussed in this section. The relative interior mass add- 
ition values of component A, f A _ B equal to 0.01 and 0.05 
are considered. For these low values of f A _ B > the result 
shows that the velocity field does not deviate much from the 
single component flows (without interior mass add- 
ition) . Moreover , the characterisitc flow features of 
driven vortices remain unaltered. Therefore, only the va- 
lues of the potential function and mass concentration are 
presented and discussed. Results of mass diffusion are 
given in such a manner that they show the effect of the 
tangential Reynolds number on the concentration distribution 
for a fixed axial Reynolds number using the pure density 
ratio as a parameter. Only the axial Reynolds number of 20 
is considered for binary flows since high axial Reynolds 
numbers are desirable to achieve convection dominated mass 
transfer which is required for fluid dynamic containment. 

It may be seen from the equation that the potential 
function is insensitive to both the pure density ratio 
and the tangential Reynolds number. This is because the 
potential function is affected by the tangential Reynolds 
number only to the extent that it affects the density 
field. And the influence of the density field on the 
potential function is mild for the relatively low values 
of f A _ B considered. This has also been shown by the results. 
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Ibr this reason, the potential function is presented only 
for a sinqlc value of the pure density ratio, 0.5, and a 
tangential Reynolds number of 100. It is given in Fig. 32 
as representative for a family of such figures. It is to 
be noted that the potential function equation is decoupled 
from the rest of the equations of the system when the pure 
density ratio is unity. In this case the potential function 
is determined once in the beginning for a specified mass 
addition schedule and the result is used for subsequent 
calculations of the remaining variables. Results given in 
this section are all for a fixed value of aspect ratio of 1. 

Results for concentrations of fissionable material 
component A, are plotted in Fig. 33, 34 and 35 for tangen- 
tial Reynolds numbers of 0, 50 and 100 respectively for the 
case of a fixed Schmidt number of 1, f A-B ®0.01 and pure 
density ratio of 0.5. It is seen in these figures that as 
the tangential Reynolds number increases the over-all level 
of concentration increases and also its gradient becomes 
steeper. The mass fraction is also plotted in Fig. 36 
through 41 for pure density ratios of 1 and 2 each for 

t 

tangential Reynolds numbers of 0, 50 and 100. These are 

also for a Schmidt number of 1 and a mass addition rate, 

f. „ of 0.01. The fiqures show little variation in concon- 
A-13 

tration distribution due to the variation in the pure 

density ratios from 0.5 to 2.0. It may be stated that 

the effect of the pure density ratio on the resulting density 


field is mild and hence the concentration equation is only 
weakly affected by the variation in the pure density ratio. 
The effect of rotation on concentration in these cases is 
similar to the case of the pure density ratio of 0.5. Both 
the over-all level and the gradient of concentration are sh- 
own to increase with increasing tangential Reynolds number. 

Next, the case with interior mass addition, f A-fi , of 
0.05 is considered. The potential function for a pure 
density ratio of 0.5 and a tangential Reynolds number of 
50 is plotted in Fig. 42. The distribution of the potential 
function obtained is similar to the case with an interior 
mass source of 0.01 since the function has been normalized by 
appropriate reference quantity. The mass fraction of com- 
ponent A is given in Fig. 43 and 44 for the tangential 
Reynolds number of 0 and 50 both for the case of Schmidt 
number of 1 and pure density ratio of 0.5. The figures 
show that the distribution is similar to that obtained with 

4 

a mass addition rate, f A _ D , of 0.01 discussed above. But 
the over-all concentration is increased by a factor approxi- 
mately equal to the ratio of the mass addition rates. It 
is also seen in this case that the vortex motion increases 
the over-all concentration as well as its gradient, which 
in. i cates the fluid dynamic containment property of driven 
vortices . 

Finally, the containment factor is calculated accord- 
ing to the following expression: 


Total Mass of Fuel (Component A) in the Chamber 
n Total Fuel Mass If 100% Fuel in the Mass Addit- 

ion Region 


P B p / / p*(r*,z*)c(r*,z*)r*dr*dz* 
* o o 

r % z$ 

p a p > I 2 r *dr*dz* 

A ' p r* z* 


( 95 ) 


The containment factor as expressed above represents a 
ratio of an actual average residence time of fissionable 
material based on steady operating conditions to that of 
an ideal reference time. This reference time is chosen 
as the residence time of fissionable material for an ideal 
case of 100% fissionable material in the assumed contain- 
ment region. The containment region here is the region in 
which the fissionable material source is uniform. 

The containment factor is computed numerically by a 
double integration using a trapezoidal rule after the steady 
state solution is obtained. The results for pure density 
ratios, p /p 0 p of 0.5, 1.0 and 2.0 are plotted against 
tangential Reynolds numbers in Fig. 45 for a mass addition 
rate, of 0.01. The figure indicates that the containment 
factor increases with increasing tangential Reynolds number 
for a fixed axial Reynolds number, Schmidt number and pure 
density ratio. The containment factor calculated for a mass 
addition rate, f A-0 , of 0.05 and pure density ratio of 0.5 
is also included in this figure. The figure shows that it 
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nearly follows the curve of f A _ 0 -O.Ol when the value is 
scaled by the ratio of the two mass addition rates. This 
indicates that the over-all concentration of fissionable 
material and hence the containment increases more or less 
linearly with the fuel mass injection rate. As may be seen 
from tiie mass diffusion equation, a further increase in 
containment factor can be realized by an increase in either 
the axial Reynolds number or Schmidt number. Finally, it 
is seen that decreasing the pure density ratio increases- 
tho containment for fixed rotation. 
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Fig. 37. Mass Fraction of Component A ( Simulated Fuel ) 
for Re = 20, Re^ * 50, Sc « 1, Pa,P^ P B,P = ^ ' 
f. - * 0.01 and a * 1. 
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Pig. 30. Mass Fraction of Component A ( Simulated Fuel ) 
for Re = 20, Re^ * 100, Sc * 1, Pa,P^B,P * 
f. « = 0.01 and a * 1. 








Mass fraction , c Mass fraction 


0.010 

0.009 

0.008 


» 0.7.5 


Axial distance, z 


0.010 

.. j i _ • j ; 

0.008 

. J_I 'I4 


; : * 

0.006 

j ' 1 ; 

' TirtT 

0.004 

. .. j. 

i • j. 


• i : i -O 

0.002 

. . ... ! i — -i._Q 

: l ' i 

0.000 

■ . i • « 

- .. . * 


•0.75 


4 0 

Radial distance, r 

Pig. 39. Mass Fraction of Component A ( Simulated Fuel ) 
for Re * 20, Re fc = 0, Sc = 1, P AfP /P B/P * 2.0, 


f* « » 0.01 and a *» 1. 
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Pig. 41. Mass Fraction of Component A ( Simulated Fuel ) 
for Re * 20# Re t * 100# Sc * X, PA,p/ p B,P “ 2 ' 
f. - = 0.01 and a * 1. 
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Fig. 43. Mass Fraction of Component A ( Simulated Fuel ) 
for Re = 20, Re t * 0, Sc « 1, P A ,p/P B ,P “ °‘ 5 ' 
f « = 0.05 and a * 1. 











Containment factor 


A,P B,P 


. ..L.j£r„ -JL .1— V- /o • - ; 

• j. . : T A * P B » p 

...... ... 1 1 — , • : • 


^A.P /p b:P “ °* 5 r ' f A -B~" P * 05 I 



A,P' "B»P 



A,P' t ’B,P 


Tangential Reynolds number, Re fc 

Pig. 45. Containment Factor vs. Tangential 

Reynolds Number for p A p/P B p * 0.5, 
1.0, 2.0 and f A _ B = 0.01, o!o5 ( The 
Case of f. « * 0.05 is Divided by 5 ) 






A130 


Pressure Calculation 

Pressure calculations are made only £or single com- 
ponent flows and they are presented in this section. As 
mentioned in Ch. IV, the pressure is calculated by the 
method of line integration for non-rotating flows. For 
rotating flows it is calculated by solving the Poisson 
equation of pressure, in which the Fourier analysis is used 
and i Courier cosine transforms are employed as an 
ef ficis omfsutational algorithm. 

The esult for an axial Reynolds number of 1 is pre- 
sented in Fig. 46 and 47 for tangential Reynolds numbers 
of 0 and 50 respectively. In the absence of rotating motion, 
the variation in the pressure is small except in the region 
near the exhaust hole where the primary flow is accelerating. 
A decrease in pressure is shown toward the chamber axis 
near the exit hole. When rotation is present, however, a 
high radial pressure gradient is formed in the bulk region 
to balance the large centrifugal force field caused by the 
angular momentum distribution induced by wall rotation. 

The uniform radial distribution of pressure shown in Fig. 

47 is in agreement with the uniform radial distribution 
of angular momentum given in Fig. 16. Near the end walls, 
the rotating motion is retarded by viscosity and hence the 
local centrifugal force is decreased. However, the pressure 
in this region is nearly equal to the imposed pressure. The 
unbalance between the pressure and the centrifugal force 
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in this region results in a large radial inward flow 
toward the axis. The mild axial pressure gradient observed 
in the end wall boundary layer and near the periphery is 
thought to be a result of the continuity constraint inter- 
acting with the momentum field. The symmetry of the pres- 
sure profile about the mid-plane, z*»0 . 5 is shown for the 
value of the ratio, Re^/Re«50. 

Results of the pressure calculation for an axial 
Reynolds number of 20 are given in Pig. 48 and 49 for 
Re t =50 and 100 respectively. These figures indicate that 
the pressure as well as its radial gradient increase as the 
tangential Reynolds number increases. The radial pressure 
gradient is seen to increase with radial distance from 
the axis. A comparison of Pig. 47 with 48 reveals a flatter 
axial distribution of pressure with lower primary flow rate 
for a fixed rate of rotation. The discussion given 
for axial Reynolds number of 1 also applies in general to 
this case of Re=20. Finally, the pressure calculation for 
an aspect ratio of 2 is shown for axial Reynolds numbers of 
1 and 20, each with a tangential Reynolds number of 15. The 
result for Re=l is given in Pig. 50 and that for Re = 20 
in Pig. 51. The pressure profile in general is similar to 
those for the aspect ratio of 1. The only noticeable dif- 
ference is the asymmetry in the distribution about the mid- 
plane. The asymmetry in the pressure profile becomes more 
pronounced as the axial Reynolds number increases for a 
fixed tangential Reynolds number. 
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CHAPTER VI 

SUMMARY AND CONCLUSIONS 


In the preceding chapters , a numerical study of dri- 
ven vortices of single component and binary fluids in a cyl- 
indrical container is described. The complete axi-symmetric 
equations of motion and mass diffusion are solved numer- 
ically by finite difference techniques. The basic differ- 
ential equations modeling the cold flow simulation of the 
advanced Light-Bulb reactor concept includes the continuity 
equation that contains a source term. Because of the 
source term in the continuity equation, the simple conven- 
tional scheme of formulating the problem in terms of the 
stream function and vorticity fails. Therefore, a method 
is developed in this study that enables the solution of this 
problem by a formulation in terms of the stream function, 
vorticity and one other function. This additional function 
is a potential function, and is introduced so that it iden- 
tically satisfies the continuity equation that has a source 
term. Another possible application of this technique is 
to unsteady compressible flow problems in which the transi- 
ent term in the continuity equation can be regarded as the 
source term. 

The governing equations are written with the unsteady 
term included in terms of the stream function, vorticity 
and the potential function. The steady state solutions are 
obtained by large time solutions of time dependent problems. 
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The. implicit numerical scheme of the ADI method is used. 

After the kinematics axe solved, the pressure is obtained ei- 
ther by lino integration or by solving the Poisson equation 
of pressure, ftourier analysis has been applied to the Poiss- 
on equation as the solution method, in which fast 
Fourier trams forms are employed as an efficient computation- 
al algorithm. An extensive study is undertaken of the 
error analysis for the numerical results obtained in this 
study. This includes the estimate of the discretization 
error and a comparison with prior studies of other investi- 
gators . 

Numerical results are obtained to show the effects 
of various parameters , in particular the axial and tangent- 
ial Reynolds numbers on the secondary flow pattern and 
the resulting containment property. Results show that in 
the absence of rotation the flow is largely through the ce- 
ntral region. When rotation is present, however, a 
high radial pressure gradient is formed in the main body 
to balance the centrifugal force. Near the stationary end 
walls, the viscosity retards the rotational motion and the 
imposed radial pressure gradient creates a large radial 
inward flow. These end wall boundary layers interacting 
with the vortex motion formed in the bulk region, induce 
recirculating secondary closed flow cells in the central 
region of the chamber. The secondary mass flow rate in 
creases with increase in tangential Reynolds number for ® 


fixed primary flow rate and aspect ratio. It appears that 
the secondary flow rate also increases as the aspect ratio 
increases for fixed axial and tangential Reynolds numbers . 

Fluid dynamic containment of the component A (sim- 
ulated fissionable material) is realized when the convection 
effect dominates the diffusion in the mass transport process 
and the mass is held within the recirculating flow cells 
by the convection currents. For a fixed axial Reynolds 
number of 20 , the results show that the secondary flow rate 
increases with the tangential Reynolds number. This en- 
hances the relative dominance of the convection effects over 
diffusion in the mass transport process. Moreover, the 
fraction of primary flow that passes through the bulk region 
of closed flow cells decreases as the tangential Reynolds 
number increases. As a result of these flow features, the 
containment factor increases as the tangential Reynolds 
number increases for a fixed axial Reynolds number, Schmidt 
number and aspect ratio. For a given secondary flow rate, 
the containment factor also increases when either the axial 
Reynolds number or the Schmidt number increases since it 
also increases the convection effect relative to the dif- 
fusion. It appears that an increase in the aspect ratio also 
increases the containment factor. Finally, the containment 
factor increases as the pure density ratio decreases for a 


fixed Schmidt number , axial and tangential Reynolds 
numbers . 

The method developed in this thesis seems to be most 
suitable and efficient for obtaining numerical solutions 
of binary flows with interior mass source. The method 
has been applied successfully to the driven vortex flow 
problems of binary fluiu, in which one of the fluid com- 
ponents is introduced by an interior mass source. The 
numerical results showed that the fluid dynamic contain- 
ment of the simulated fuel component/ introduced as an 
interior mass source/ can be achieved by the convection 
c u rrents of recirculating closed flow cells formed by the 
rotating periphery. The containment property is shown to 
increase with increasing the rotation for fixed primary 


flow rate. 
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ABSTRACT 


The gas core nuclear reactor concept was Introduced 
fifteen years ago because of a possible application as an 
engine for space vehicles. Recently, interest has grown in 
applying this concept for ground based power generation. 

In this thesis the results of a parametric study on the 
entrance flow region in a gas core nuclear reactor are pre- 
sented. The physical system is modeled as laminar confined, 
coaxial flow with heat generation in the inner fluid. The 
governing equations include the boundary layer approximations 
and the assumptions of only radial radiative transport of 
energy represented as an energy diffusion term. 

The Von Mises transformation and a c transformation 
are used to transform the equations into nonlinear nonhomo- 
geneous convective-diffusion equations. A unique combination 
of forward and backward difference equations which yields 
accurate results at moderate computational times, is used 
in the numerical method. 

Results show that the rapidly accelerating, heat generat- 
ing inner stream actually shrinks in radius as it expands 
axially. This conclusion is in opposition to that assumed 
in previous analysts of the flow region downstream of the 
entrance region. 
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CHAPTER I 
INTRODUCTION 


Fluid mechanics problems associated with the concept of 
a gas core nuclear reactor have been studied for many years. 
Most of the theoretical work has been limited to the main 
portion of the cavity where severe gradients of velocity, 
temperature, and pressure do not exist. The reasons for 
this limitation are obvious. There is also a general lack 
of understanding about the behavior of the transport mech- 
anisms in entrance flows. Because of these problems-^ the 
highly complex entrance flow problem of gas core reactors 
has been essentially ignored. Only one previous attempt has 
been made txr solve the complex entrance flow problem and 
that analysis used inviscid flow equations. 

The object of the present work effort is to formulate 
a realistic model of the entrance region flow in the reactor 
and to study the effects of various parameters on the flow 
field. A new transformation is developed to solve not only 
this problem but entrance flow problems in general, since 
simple solution methods to this problem are not generally 
available. A unique combination of forward and backward 
difference equations is used in the numerical method, which 
yields accurate results at moderate computational times. 

In an open cycle gas core nuclear reactor concept, it 
is presently proposed that a solid uranium rod be continu- 
ously fed into the reactor cavity to produce makeup fuel for 
that lost through the nozzle and to maintain a critical mass. 

BIO 
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The fuel rod enters the cavity on the axis of the working 
fluid flow. It Is exposed to the high neutron flux and Is 
very rapidly vaporized. The temperature of the uranium 
vapor Increases by an order of magnitude In a very short 
distance downstream of the entrance point. The Interaction 
of the rapidly expanding uranium stream with the surrounding 
working fluid stream Is difficult to predict. It has an 
Important Influence on the mixing of the two streams, as 
well as on the critical mass requirements In the cavity, 
since the fuel stream must expand to a radius many times 
that of Its Inlet radius to keep critical mass requirements 
reasonable. It Is also the object of this work to analyze 
the flow In the region of the fuel Injection, point where 
temperature changes, and hence velocity and density changes, 
are extremely large. This study will provide Input for an 
analysis of the main portion of the reactor cavity where 
gradients of the flow variables are less severe. 

The analysis presented here Is an approximate one and 
attempts to consider only those phenomena which are of funda- 
mental importance in the region of Interest. The flow Is 
considered to be coaxial, with the two streams Immiscible. 
Thus, the diffusion equation can be discarded, and the prob- 
lem can be treated as a two region problem. Radial body 
forces are neglected. Axial body forces can be Incorporated 
in the pressure term. Thermal diffusion terms, except for 
Rosseland’s approximation for radiative energy transport, 
are neglected. Internal heat generation occurs In the fuel 
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stream to represent the fission process. A single equation 
of state Is used for the whole region. Initial velocities 
are very small , and velocities throughout the cavity are 
very low subsonic. Physical properties are assumed to be 
uniform. 




CHAPTER II 


BACKGROUND 

/ 

This thesis is concerned with the analysis of flow in 
the entrance region of a gas core nuclear reactor. Previous 
work relevant to this thesis work can be separated into two 
categories. The first of these is the literature on model- 
ing of flows in reactors and the second is the literature 
concerned with entrance region flows in general. 

Reactor Modeling 

Ever since the concept of gas core reactors for nuclear 
rocket engines was proposed in 1953, its development has 
been the basis for several experimental and theoretical in- 
vestigations. The feasibility of a gas core reactor has 
been studied ir» detail 12 * 1 ^ A status of the art report 
by Rom 1 ^ refers to the basic work completed in the fields 
of hydrodynamics, gaseous radiant heat transfer, neutronics 
and system studies. 

The exchange of energy from the fissioning fuel to the 
propellant stream by thermal radiation has been treated by 
Kascak^® using a diffusion (Rosseland’s) approximation to 
the radiant heat flux, which is valid for optically thick 
fluids . This treatment is adopted for the present analysis. 

The analytical work on the fluid mechanics problems of 
gas core reactor has been limited to the main portion of the 
reactor cavity where gradients are not extreme. These 
treatments were numerical integration of the boundary layer 
equations, or full Navier Stokes equations for coaxial flow 
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of dissimilar fluid3. They were concerned mainly with the 

mixing of the two fluid streams and hence concentrated on 

♦ 

the diffusive terms rather than on the effects of the large 
Internal heat generation. References 6, 7, 18 and 24 are 
examples of this approach. 

Ohia^ et al. studied isothermal, confined laminar mix- 
ing of dissimilar (binary and nonreacting) circular axisym- 
metrlc Jets. Their results clearly indicate that the dif- 
fusion is negligible for at least one characteristic length. 

pp 

Only one previous attempt** was made to simulate the 
flow field near the fuel injection location. Using inviscid 
equations this analysis showed the effects of large energy 
release and the concurrent high acceleration of inner stream 
fluid on the coaxial flow field. Results of this study in- 
dicated that the rapidly accelerating, heat generating, inner 
stream actually shrinks in radius as it expands axially. 

Even though inviscid analysis yielded certain useful results, 
it can not simulate a flow with interfacial shear and the 
effects of solid wall on the flow field realistically due to 
the degenerate nature of the equations. Also, in the invis- 
cid analysis, stability problems were experienced and 
artificial viscosity type terms were added to the governing 
equations to overcome these difficulties. Based on this 
previous experience the present analysis is performed using 
the boundary layer equativ .3. 


general Entrance Plows 


There has been an enormous amount of published liter- 
ature on a variety of entrance flow studies. Unfortunately, 
none of these is directly applicable to the present analysis. 
A general description of this literature Is, however, given 
to Justify the methodology used in formulating and solving 
the present problem. It is impractical to cite all the pub- 
lished work on entrance flows; therefore, only recent papers 
dealing with circular conduits which contain new ideas 
and/or a review of previous work are cited. 

The reason entrance region flows have been studied so 
intensively is because of their practical importance in the 
investigation of excess pressure drop, flow stability, and 
flow separation phenomena. Most of the studies involve 
solution of simplified Navier Stokes equations using bound- 
ary layer approximations. However, there has been some 
controversy over neglecting the axial diffusion terms in the 
governing equations for small Van Dyke showed rigor- 

ously in the case of entry flow in a channel that the full 
Navier Stokes equations reduce to boundary layer equations 
for the infinite Reynolds number flow. (Solutions of the 
full Navier Stokes equations^ indicate that the boundary 
layer approximations are valid only for Reynolds numbers 
greater than 250.) In conclusion the flow field description 
obtained by using boundary layer equations is reasonably 
accurate for the asymptote of large distance from the 
entrance when Reynolds number is large. Boundary layer 
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equations are preferred for the present analysis because 
general flow features are the main concern. Also, the 
methods of solution for the full Navler Stokes equations 
become cumbersome and/or have stability problems for large 
Reynolds number flows. 

Pour basic methods of solution for the boundary layer 
equations are encountered In the literature: 

1. Finite differences of the complete set of 
equations 

2. Linearization of Inertial terms 11 * 1 ^; 

3. Method of weighted residuals 1 ; 

li I? 

and 4. Patching of solutions * 

Most of these methods are numerical and all have their own 
merits and drawbacks. Even though analytical solutions are 
preferred, numerical techniques are widely used in engineer- 
ing practice because they have the capability of handling 
more comprehensive cases with variable fluid properties and 
variety of boundary conditions. The simpler numerical 
schemes, even though they are easy to set up, usually become 
unattractive because large numbers of grid points are required 
to resolve the flow characteristics near solid boundaries and 
in regions of large gradients. Also, the simpler methods 
have stability and convergence problems when derivatives are 
large. This problem can sometimes be overcome, however, by 
judiciously introducing transformations Into the basic set of 
equations, which results in new independent and dependent 
variables before using the numerical technique. Obviously, 


this procedure will work when derivatives of the new depend- 
ent variables are much smaller than those of the original 
dependent variables. 

In this thesis a series of transformations are used 
before solving the partial differential e quatio ns , A 
combination of explicit and implicit finite difference 
methods is also utilized with the transformed variables to 
obtain the largest possible grid size consistent with 
stability and convergence criteria. 


CHAPTER III 


MATHEMATICAL MODEL 

The geometric model for the analysis consists of two 
fluids flowing coaxially within a bounding stream tube of 
constant radius (See Pig 1). Radial body forces are neglect 
ed to preserve axial symmetry in the flow system. Under 
these assumptions the steady flow situation is represented 
mathematically by boundary layer equations with appropriate 
boundary conditions. The inner stream has a large internal 
heat generation rate and transfers heat to the outer stream 
by radiation. The radial conductive transport of energy is 
negligible com p a red -to - -the--radl a l-r adlatlve transport. A 
diffusion approximation is used for the radiation equations 2 
The fluids are considered to be immiscible and no mass trans 
fer occurs across the interface between the two streams. 

The interface is a stream line of both fluids. Since the 
velocities, temperature and pressure keep changing as the 
fluids travel through the tube, the radius of the interface 
does not remain constant. The equations are written for 
both fluids separately and tied together with the interface 
equation. 

The following are the governing differential equations: 
Fluid I: 0<r<_aR. j = l 

Fluid II: a R < r < R. .1 * P 

State: Fj » itj pj Tj ( 1 ) 

§? (f J 5 J ?> * I? <5 J Vj F) = o 


Continuity: 


( 2 ) 


CONFINING CONDUIT 



Figure 1. System for Coaxial Entrance Flow Study 





Momentum: 







•a* 

a? 



+ 



a 

3r 


_ 3u* 

(r -=?-) 
3r 


( 3 ) 


Energy : 


_ _ 3T< 3T 1 x 

>, c 1 (Ui + v 1 - 3 *-) 

j j >1 j z J 3r 


Ki ft 3 3T« 

5 j s j *i + i h (F t 3 #> 


Here the inner stream is assumed to generate heat at a 
constant rate per unit mass. For the case when the inner 
stream generates heat at a constant rate per unit volume, 
equation (4) can be written as: 


_ _ 3Ti _ 3T 1% A ki a _ _3 3Tj 
>4 c, (u. — + v. ~zr) * 4*4 + — ( r T. —zr~) 

J J J 3 Z J 3r J r 3r J 3r 


From now on the subscript j, indicating either the inner 
fluid or the outer fluid, will be omitted for simplicity. 

The following are the boundary conditions at the wall 
of the tube: 
u * 0 


v * 0 and (5) 

T ® T . 

* *w* 

Due to the axial symmetry of the problem, the condi- 


tions at the center line can be written as: 



(aR) * U 2 (aR) 


du-. 

Mi — 

1 3r 


k l T 1 — 

1 dr 


r * aR 
3 3^i| 


3u- 

- u 2 -=r 
^ 3r 


-3 3^2 
k 2 t 2 — 


aR 


( 8 ) 


| r * aR I r * aR 

The initial conditions are specific to a given calculation. 

To nondimensionalize the equations, the following 
dimensionless variables and groups are introduced: 


r 



’M 



z 


z 

R Re 


Re 


V J 


2i 


-? 


•> • « 



4 — 


. 1 ? 



" *S 


u & 


T. * - 
J T 


r J »c 
P 2i u ?l 


R u 2i P 2 i 


k 2 t 21 
©2 U 2 


*1 ■ — Re Prl/3 

J c 2 I 21 


*.1 6 c T 2i 
4^173 


Substituting the above into the dimensional equations 


we obtain: 


Fluid I: 
Fluid IT: 


State : 


0 < r < o, 
o < r <_ 1, 
P » 0 p T 


j * 1 
J * 2 


3 3 

Continuity: ■yj (p u r) (p v r) a 0 


( 10 ) 


Momentum: 


„ „ 3u , . 3u 
p u 3z * p v 3r 


+ M ?f? 


<'$?> 


( 11 ) 


Energy: 


3T 3T 

PU-g^+PVjp«PU* + 


§ F !r T? !?> 


( 12 ) 


M 
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The boundary conditions are: 
u (1) ■ 0 
v (1) - 0 

T (1) - T e 
3u(0) 


Tv 
v (0) 


0 

0 


3T(0) _ n 
3r “ u 

and the interface conditions are: 


(13) 


(lH) 


da _ V 1 
dz " u^ 

r « 

v 2 

m — 

U 2 

a 

r * a 

(15) 

P x (a) 


- p 2 

(a) 


T x (a) 


* T 2 

(a) 


u 1 (a) 


* u 2 

, (a) 

(16) 

M, 

1 dr 

3u 2 (a) 

3r 


K -i (a> 

K t 3r 

3T 2(a) 

3r 



The equation set (9)-(12) and (15), with boundary conditions 
(13) and (1*0, and interface conditions (16), is well defined 
with a proper set. of initial conditions. 

The Von Mises transformation relationships 

—£• ** pur 
3r 



are Introduced to reduce the equations In number by elimi- 
nating the continuity equation and the radial velocity 
component, v, to the following equations in (z,V>) coordinates 


pu lz + I§* Mpu l^ (pur 2 |^) 

3T K 3 , 2 _3 3T. 

Tz m S 7* (pur T 


( 18 ) 


(19) 


The equation of state Is unchanged under the transformation. 
It is Interesting to note the equation (15) represents a line 
In the (z,ty) plane. 

Though the equations become simpler In Protean coordi- 
nates, three serious difficulties come into focus; 
i) The boundary condition at the center line in (z,<|>) 
plane Is 


3u 

1 ? 


• 1 3u 0 

* pUF 5r * 


|r * 0 

Using L'Hospltal's rule this may be evaluated as 


3u 

pu ?* 


Similarly, |X| 
3t|» 


3u 

3r* 


* finite 

r » 0 

can also be shown to be bounded 


and nonzero. Notice that symmetry Is lost with respect 
to ip at ■ 0. 

11) The first and higher derivatives of u with respect to 
<j» at the wall are Infinite. This means that wall is a 
singular line. Hence Taylor Series expansions cannot 
be written at the wall. Moreover, the truncation 
errors in finite difference molecules will be rather 
large near the vicinity of the wall, 
ill) If a uniform finite difference grid is used in the 

plane, this would amount to having a fine grid 
near the center line and a relatively coarse grid near 
the wall in the (z,r) plane. Hence, a large amount of 
points are needed to simulate the flow situation with 
reasonable accuracy. 

To overcome these difficulties, the following trans- 
formations are introduced: 


W 

0 

♦ 

L 


U 

2 


T_ 

4 

.a 


2 


1 4 

- - ; ) l 

u l i P 11 r dr * 'J u 2i p 2i rdr 


2 /j- (1-r 2 ) rdr 


( 20 ) 

( 21 ) 

( 22 ) 


L Is the ratio of total mass of both fluids to the muss 
of an incompressible fluid entering with a uniform dimension- 
less inlet velocity of 1 and a dimensionless inlet density 


M “ 
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of l. The use of a scale factor, L, Is convenient, when 
simulating coaxial flow of two Incompressible fluids. The 
methodology used In arriving at the c transformation Is out- 
lined in Appendix A. 

Equations (18) and (19) transform Into the following 


equations with some manipulation: 

2 

aw . i dP _ i / v 3w . & awt 

5¥ + pdi°E {Y ^ +6 I? } 


(23) 


30 


- H 


V 


{a 


d 2 e 

a? 


30, 

+ “H> 


(2*1) 


where 


JS 


u fC 
h u 


JL £1 


u 4 


|X 


H = 


* t 3 


S Uf 
K T 3 ' 


a = y u 


u> = u 6 + I 

U <K 


and 


u 


= 2(1-C 2 ) L 


(25) 


v 




T 


It should be noted that equations (23) and (2<l) are 
analogous to a pair of coupled unsteady nonlinear nonhomo- 
geneous diffusion equations. Also, u^ Is the solution to 
the equations describing fully developed flow of two coaxial 
streams of Incompressible fluids with no generation. 

The scaling factor, L, accounts for the effect of 
different entrance velocities and densities of the streams 
on the magnitude of the fully developed velocities. The 
use of L results in mapping of the region 0 < r <_ 1 into 
0 < 5 < 1 , 

The boundary conditions for the above set of equations 


transform into: 
W(l) - 0 


( 26 ) 


0 ( 1 ) 


inf + 


(27) 


3W(0) 


30(0) 

3 ? 


( 28 ) 


(3) - P 2 <*> 
© x (3) » © 2 (3) 
W x (3) - W 2 (3) 


3Wi (g) Q X 3W 2 (0) 


‘1 FT 


o 2 H 


O'* 


o 
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K !fl (0) . jjl 80 2 (3) 

1 3? ^2 H 

where 8 is evaluated from the following equation: 


( 29 ) 


r o u f « 


'o Pj Uj r dr 


(30) 


Notice that symmetry at the center is preserved under 

the transformation and c 8 1 is not a singular line in (z,e) 

plane. The condition on P [equation (27)1 is arrived at by 

using the momentum equation at the wall with the boundary 

condition (13) and applying L’ Hospital’s rule to evaluate 
3u 

the limit of at the wall, 
r 

Obtaining an analytical solution to the above system of 
equations is not possible due to the highly nonlinear nature 
of the problem. Hence, the solution is obtained using 
numerical methods. 

The fully developed solution for the entrance flow 

problem is obtained by solving the equations resulting from 

3u 3T 

dropping the terms and from equation set (9) through 
(16). The fully developed velocity profile is given by 

„ . a 2 + ”l (l-i » 2 ) - r 2 

dz M, 

0 < r < a 



a < r < 1 


The fully developed temperature profile cannot be 
obtained analytically for the case where the inner fluid 




generates heat at constant rate per unit mass; but for the 
case where the Inner fluid generates heat at a constant rate 
per unit volume It Is given by 
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T e * + *S t(a 2 -r 2 ) - ^ ln o) 


0 < r < a 



2 

♦S SL. In r 
2 


a < r < 1 


In either case the values of ~ and radius of inner stream 
for fully developed flow cannot be obtained analytically. 
However, ln the Incompressible Isothermal case with two 
fluids of the same viscosity the following relationships can 
be derived: 


dP 

dz 


and 0 * 


- 8 (u <*2 ♦ 1-or ) 
R 

/ 2ul 

>‘V‘- -V s 


u R o 2 + 1-a 2 


When pp » 1 and Up a 1 these relationships agree with 
those for the classical entrance flow through a circular 
conduit. 
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CHAPTER IV 


NUMERICAL METHOD 

Several numerical methods are available for the solution 
of a parabolic system of coupled nonlinear partial differen- 
tial equations. The finite difference methods were used 
because of their simplicity. Also, the stability, conver- 
gence and consistency for these methods have been very well 
established. In the finite difference technique two proce- 
dures exist, the explicit, and the implicit. In the explicit 
methods, forward differences are used to approximate the 
derivatives in the axial direction, which lead to linear, 
explicit finite difference equations that can each be solved 
independently. The implicit methods, however, use backward 
differences to approximate the axial derivatives, resulting 
in a system of simultaneous, nonlinear equations. In 
general, the implicit methods are unconditionally stable, 
whereas the simple explicit methods are inherently unstable. 

In order to obtain stable, accurate solutions by using the 
explicit methods, restrictions on the step sires are necessary. 
The proper choice of a computational method depends on the 
nature of the problem at hand. 

The application of a purely explicit or a purely 
implicit method to the system of equations describing the 
entrance flow was not satisfactory. Due to the nonlinear 
nature of the problem, the restriction on maximum allowable 
step size to ensure stability in the case of an explicit, 
method became progressively worse and the computational time 


B31 


became prohibitive as the temperatures of the fluids started 
increasing. On the other hand, convergence problems were 
experienced in the solutions of simultanr us nonllrear equa- 
tions near the entrance region in the case of implicit 
methods. Hence it was decided to use the explicit method 
near the entrance region and change to the implicit method 
after the gradients become less severe. This approach yields 
highly accurate results over the whole region at moderate 
computational times. 

The flow region is descretized by grid points having 
equal spacings (o<c<0) and Ac 2 (8<;<1) in the c direction 

and an arbitrary axial step size Az. The discretized rect- 
angular grid and the coordinate system used to solve the 
problem are shown in Figure 2. The subscript n and the 
superscript i are grid point indices associated with C and z 
respectively. The radial subscript n takes integer values 

between 1 and (M+N+2) and the axial subscript i ranges from 

♦ 

1 = 0 to any desired value. The interface is labeled by 
u = N+l for the inner fluid and n * N+2 for the outer fluid. 

PI nlt e Pi fference Equations - Explicit Method 

Substitution of three point central difference approx- 
imations for the radial derivatives and of two point forward 
difference approximations for the first derivatives in the 
axial direction into Equations (23) and (2*0, results in the 
following two finite difference equations for every interior 
point, (I, n) of the finite difference network: 



Finite Difference 3ri3 System 



1 


W n- 


4P i Az r .A n+1 

—t + — r It* 

P 1 E 1 n 


- 2W 1 + W 1 , 
n n- 1 


4 . At n+1 - W £-l , 

n TnI > 


(31) 


3* + H Ax + {o* 0 ^>±l „-^ e n + ° n -l 


JL i 

+ u i e n+l ~ e n-l , 
n 777 ) 


(32) 


The subscript j , denoting the first or the second fluid, 
has been omitted for simplicity. 

The numerical stability considerations impose the fol- 
lowing restriction on the largest value of Az that can be 
used in the above scheme for the known Ac 111 . 

Az < Min { | I A; 2 1 , I ° a; 2 | ) (33) 

lall n c r lall n 

The value of Ac is not limited from stability con- 
siderations and is selected from the required resolution 
and the accuracy of the flow problem. 

The boundary conditions, given by equations (?b) 
through (29), have the following finite difference 
representations : 


W (M+N+2) 
0 (M+N+2) 


* 0 \ 

- 0_ / 


(34) 


3U 3u , 

AP 2 " p M+N+2 *2 { — 2 + r-l) 

’«• *«• In-l 


(35) 


If 



derivatives at the center line, the interface, and the wall, 
one sided differences of second order accuracy have been 
used. 

After prescribing a set of initial conditions, the 
following sequence of operations is performed to obtain the 
solution in the entire region of interest : 

1. Compute Az for all n from the stability condition 
(33). Use the smallest value of all these values for further 
computations . 


1 


■■'•-T7T 


““1? - Vjv 


;V w 

V ,/'b 
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2. Compute the value of AP from equation (35). 

3. Compute the new values of W and 0 l’ron equatlonn 
(31), (32), (33). (35) and (36). 

4. Evaluate tho location of stream lines and values 
of radial velocities by solving the following equations: 



(38) 


3r\ 


(39) 


5. Proceed to next axial step. Repeat steps 1, 2, 3 
and 4. 

Finite Difference Equations - Implicit Method 

The finite difference equations resulting fr»r>m the 
substitution of central difference approximations of fourth 
order accuracy for c derivatives and of two point backward 
difference approximations for the first derivatives in the 
2 direction are: 

(16RY* +l -8T 1+1 ) W i+ J - (30RY i+1 +E i+1 ) W i+1 
n n n-i nn n 


4* (16RYf +1 +8T6^ +1 ) W 1 !} 
n n n+1 


E l+1 


n 




• w i> 


( R Yn +1 


* T«J +1 ) 
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0 


0 


(16RaJ + l-8Tu»J +1 ) ©*+* - (30Ro*+l + F* +1 ) ©J+ 1 

+ (16Roi +1 +8Twi +1 ) G 1 *} 
n n n+i 

(41) 

. - (©i+H^Az) F 1+1 + (Ro i+1 -Tw i+1 ) 0 i+ i 
n n n n n n- 2 

+ (Ro* +1 +Tu* +1 ) ©ft* 
n n ntz 


0 


R 


Az 

12 (AC ) 2 


(42) 


m 


0 


T 


Az 

12(4C ) 


(43) 


The finite difference representations for the boundary 
conditions are given by equations (34) through (37). Equa- 
tions ( 36 ) and (37) take a slightly different form because 
of higher order approximations. For any grid size, the use 
of five point finite difference approximations for deriva- 
tives results in improved accuracy, compared to three point 
finite difference approximations. 

The following procedure is used to advance the solution 
from the i th row (where W and 0 distributions are known) 
to the i+l th row: 

1. Assume a value for AP. 

2. Assume (N+M) values of W and (N+M) values of G for 
the i+1^* 1 row. 


.? 
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ft 
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3. Compute new values of W and 0 by solving N+M equa- 
tions of the form of Equation ( 40 ) and of Equation (41) by 
tridiagonal matrix Inversion. 

4. Use the newly generated values of W and 0 and 
modified form of Equation (36) to evaluate new values for 
wj+l and 0* +1 ; also compute new values for location of stream 

lines from Equation (38) and a new value for AP from Equation : 

(35). 

5. Test AP, W and 0 for convergence of the Iteration 
scheme. 

6. If the convergence requirements are not met, repeat 
steps 3, 4 and 5. If the convergence criteria are satisfied, 
compute the values of v and proceed to the next axial step 

in the marching procedure. 

i 

With this computational scheme the solution is initi- ' 

ated at the tube entrance, z * 0, where the flow field Is j 

prescribed through the initial conditions, and is then obtain- j 

od at each successive axial position marching down the tube. j 

, . . I 

Though there is no restriction on the maximum value of Az j 

used in the computation from numerical stability considera- 

/ 

♦Ions, the set of initial conditions prescribed and the 
accuracy required in the final values dictate the practical 
value to be used. In step 4 of the Implicit method computa- 

i 

tional scheme, to obtain an improved estimate of the value 

J 

of AF, the use of * ~laxation and gradient techniques was j 

proven to accelerate the convergence of the iterative scheme. 

♦ 

\ 

l 





CHAPTER V 
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RESULTS AND DISCUSSION 

The results of a parametric study of the coaxial en- 

* trance flow with heat generation are Illustrated In Figures 
3 through 2?. The values of various system parameters 
choosen for run No. 1 are listed In Table 1. For this run, 

* heat generation was prescribed as constant per unit mass. 
Initial velocity and temperature profiles were assumed to 
be flat. The base run has an axial velocity ratio of inner 

* to outer streams of 0.1, and a density ratio of 10. Both 
streams have the same initial temperature. The values of 
specific heat, thermal conductivity, and gas constant for 

* uranium and hydrogen were used in evaluating other dimen- 
sionless groups. 


Table 1. System Parameters 



Parameter 

Value 

0 

4 

6. 333 


k 

0.01 


S 

0.01 


°i 

1.776 (107) 



1.776 (10 8 ) 


a 

0.7 

* 

Pr 

1.0 


T 

e 

1.0 


UJ " - ( jr* 





* 9 


SITunit iff* 1 /* 10 ® 1 ?' Profllea ~ Constant Oenerat 
par unit Nmi Cm* - Paraneter Set Vo. l 





Flexure 4. Temperature Profiles - Constant Generation 


per Unit Mass Case - Parameter Set No. 1 

















Figure 7. Temperature Profile# - Constant Generation 
per Unit Mass Case - Parameter Set No. 2 










Figure 9. Axial Velocity Softies - Constant Generation 
per Unit Maas case - Parameter Set No. 3 

















figure lfc. Axial Developnent of Central and Interfacial Axial Velocities and 
Temperatures , Pressure Drop and Location of Interface - Constant Generation 
per Unit Mass Case - Parameter Set Bo. k 




Figure 15. Axiel Velocity Profiles - Constant Generation 
per Unit Maas Case - Parameter set No. 5 



Figure 16 . Temperature Profiles - Constant Generation 
per Unit Mass Case - Parameter Set No. *? 











Figure 1 M . Temperature Profiles - constant Generation 
per Unit Volume Case 
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A total of seven runs were made to study the effect of 
various parameters on the flow field. For the first five 
runs, the generation term was chosen as constant per unit 
mass of fluid. The parameter .riations for the first five 
runs are listed in Table 2. The sixth and seventh runs were 
made to study the effects of different forms of generation 
terms. For the sixth run generation term was chosen as con- 
stant per unit volume. Coaxial entrance flow without heat 
generation was simulated by the seventh run. 

Table 2. Parameters for Constant Generation Per Unit Mass Case 


Set No. 

Parameter 
which is 
Varied 

Value 

1 

Base 

- 

2 

U R 

0.05 

3 

PR 

20 

4 

k 

0.005 

5 

a 

0.5 


Before discussing the results obtained, it is necessary 
to explain the shape of the initial profiles chosen. 

Initial Profiles 

Flat initial profiles were chosen for simplicity in 
calculation. Even though a singularity exists where the 
Inner stream meets the outer stream, this flow would be more 
stable than two parabolic profiles meeting at the inlet of 
the reactor. 







V. 
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There has been published experimental data which indi- 
cate that a back flow region Is formed at an outer to inner 

stream velocity ratio of 2(i for heterogeneous Jets with a 

X 6 

density ratio of 4 . Since the nature of the parabolic 

governing equations is r.uch that they cannot handle flow 
situations where there is a substantial radial variation of 
pressure which will eventually lead to the axial or primary 
Plow velocity component becoming negative, a value of 10 is 
chosen for the outer to inner stream velocity ratio. This 
ensures physical stability in the computational scheme, 
since it limits the size of the radial pressure gradient 
to a negligible fraction of the axial pressure gradient. 
Validation of Results 

There are no known experimental data for the laminar 
confined coaxial entrance flow problem with or without heat 
generation. Therefore, checking the validity of the re- 
sults by direct comparison is not possible. However, 
verification was attempted by simulating the classical iso- 
t. hernial , laminar, Newtonian flow in a tube entrance region 
for a uniform entrance velocity. Computed velocity profiles, 
entrance lengths and pressure gradients are in excellent 
agreement with those reported in published literature. 

Details of this comparison are described in Appendix B. 
l.ven though this proof of validity is not conclusive since 
it Is inferred from simulating a special case, it provides 
a partial check of the method used. 
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Discussion of Results 

The developing axial velocity profiles and temperature 
profiles are shown In Figures 3, 4, 6, 7, 9, 10, 12, 13. 15. 

16, 18, 19 and 21. Figures 5, 8, 11, 19, 17, 20 and 22 show 
the axial development of center line and Interfacial axial 
velocity and temperature, Interface location and pressure 
drop. All the results reported were generated using a 22 
point c grid with 11 points In each fluid. The explicit 
method finite difference scheme was used up to s - 0.03 and 
the implicit method was used thereon until s * 1.0. The 
convergence of the solution of finite difference equations 
to the solution of partial differential equations was proved 
by comparing the results obtained for run No. 1 of constant 
generation per unit mass case using a 42 point t grid, with 
those obtained using a 22 point C grid. The results match 
within 1 % accuracy. 

The axial velocity profile plots Indicate that in all 
cases the flow exhibits the same type of behavior. Close to 
the entrance, the faster moving outer stream transfers some 
„f Its momentum to the Inner stream, thus accelerating It. 
Also, the effect of the outer stream coming in contact with 
the wall, which forces the velocity of fluid at the wall to 
be zero, Is manifested by a peak In the outer stream axial 
velocity to satisfy continuity. The maximum observed in the 
uxtal velocity keeps traveling towards the center line as the 
flow develops axially, and In all cases It Is seen that at 
2 = 1.0 the velocity profiles are parabolic. 



The results show that the Inner stream shrinks as it 
travels through the tube. This can be explained singly 1" 
the case of sere heat generation. The outer stream con- 
tinuously transfers some of its momentum to the inner stream, 
thus accelerating It. The inner stream must contract in 
radius to satisfy continuity. This happens until the flow 
is fully developed. 

The internal generation of heat in the inner stream 
raises the fluid temperature. Since the specific volume of 
the fluid is proportional to the temperature, the value of 
the center line velocity is further enhanced. This again 
could result in a shrlnking-while-expanding inner stream. 

A reasonably accurate prediction as to whether the inner 
stream is going to expand or shrink in radius can be made 
from the value of the ratio of the product of average density 
and average axial velocity at any z to that at z - 0 of the 
Inner stream. The inner stream will shrink or expand if 
the ratio is greater or less than unity, respectively. All 
the runs reported showed a shrinking inner stream. Attempts 
to obtain values of parameters necessary to generate an 
expanding inner stream were abandoned due to severe numerical 
II faculties. It is felt that full Navier-Stokes equations 
are necessary to describe such a flow situation. 

The profiles generated and flow features are in good 
qualitative agreement with those reported using lnvlseid 

2 ? 

flow analysis . 
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KftVcts of Plow Parameters 

PiKur»'n 5, 8, 11, 14, 17, 20 and ?i present the offects 
of the following flow parameters on flow variables: 
l. Magnitude and type of generation term, 

<■> u ratio of inner to outer stream initial entrance 

velocities, 

3. p R ratio of inner to outer stream initial entrance 

densities , 

4. K ratio of inner to outer stream thermal con- 
ductivities, and 

a, dimensionless Inner stream radius. 

Figures 5, 20 and 23 show the effects of magnitude and 
type of generation term on the flow variables. With no 
generation, values of center line velocity u x and inter- 
racial axial velocity u N+1 reach values of 1.11 and 1.06 
respectively, at z * 1.0. The inner stream has a diameter 
a of .21 and the pressure drop DP has a value of 5* 47. With 
the same flow parameter values, but different generation 
terms, there is a significant increase in the values of all 
variables Just described. In the first case where the inner 
stream generates heat at a constant rate per unit mass, 
u u N+l» cenler llne tera P erature T i» mterfacial temperature 

T a and DP attain values of 1.84, 1.66, 4.39, l»8l, .31 
and 9.04, respectively, at z « 1.0, compared to 1.4, 1.28, 
2.93, 1*3* *29 and 6*78 respectively, in the second case 
where the fluid is assumed to generate heat at a constant 

Values in the second case are lower 


rate per unit volume. 
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than those in the first case due to the smaller amount of 
heat generated. If the inner stream is expanding it la 
expected to result in an Increase in all values. A com- 
parison of values of product of average density and average 
axial velocity explains why the diameter of the inner stream 
is larger in the first case compared to the second case. It 
is also interesting to note that the inner fluid temperatures 
reach a maximum near the entrance for the second case and 
gradually become smaller as the flow develops. This is 
completely opposite to the behavior in the first case. A 
close look at the phenomenon explains this behavior. 

Initially, the inner stream occupies a large area and as it 
goes down the conduit it starts to shrink in size and hence 
the amount of heat generated decreases. Since the fluid 
cannot generate enough heat to compensate for heat transferred 
1.0 the outer fluid, thereby sustaining the initial tempera- 
tures, the temperatures start falling off. Similar behavior 
is expected to be seen in the first case also, if the magnitude 
nf generation term is not as large. 

The effect of change of other parameters is only report- 
ed for the constant generation per unit mass case. A similar 
study was also made for the constant generation per unit 
volume case, but the results are not reported since no new 
Information was obtained. 

By decreasing the value of u R to 0.05, the final values 
of u N+1 , T x , T n+1 , a and DP decreased to 1.54, 1.47, 3.7, 
1.7, .22 and 6.98, respectively. This is due to the decrease 


i 
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In total Initial momentum and amount of heat generating 
fluid. The opposite reason explains the lnorease In the 

values of Uj, U N+1» T l» T N+1» ° and DP t0 2,1 » 5.22, 

2.12, .32 and 12.15 resulting from an increase In the value , 

of p R to 20. A decrease in the value of K reduces the amount 
of heat transferred from inner stream to outer stream, 
resulting in a hotter inner stream and a colder outer stream. 

This leads to smaller u^, u N+1 and DP values, and a larger 
value of a. The temperature profiles exhibit steeper 
gradients in this case. A reduction in the amount of inner 
fluid entering the tube (i.e., smaller value of a) results 
in a faster moving inner stream (due to increased total 
initial momentum) with smaller temperatures (due to lesser 
heat generation). 

A comparison with the fully developed profiles was not 
made since the flow does not seem to have been fully developed 
at z 8 1.0. 

i 


> 


1 
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CHAPTER VI 

CONCLUSIONS AND RECOMMENDATIONS 

The original contributions of this theala are a model 
describing transport mechanisms In the entrance region of 
the gas-core nuclear reactor, and a c transformation which 
makes possible the use of a protean coordinate system for 

confined flows. 

The results obtained In this thesis using the laminar 
flow model show that the rapidly accelerating, heat generat- 
ing inner stream actually shrinks in radius as it expands 
axially* This conclusion is in opposition to that assumed 
in previous analyses of the flow region downstream of the 
entrance region. 

Since the attempts to obtain an expanding inner stream 
by varying the basic system parameters were not successful, 
it is recommended that full Navier-Stokes equations be solved 
to gain better understanding of the complex entrance flow 
problem. 

The value of c transformation as an efficient method to 
solve confined flows should be further explored. 


APPENDIX A 
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C TRANSFORMATION 

Finite difference schemes using grids with uniform 
spacing are the simplest and the most convenient for program- 
ming. However, such schemes are not satisfactory for use in 
problems with boundary layers. The numerical solution could 
have gross errors even in the stream core, if the number of 
points is riot great enough to resolve the boundary layer. 

The use of enough grid points to resolve the boundary layer 
makes the computational time unacceptably large. This prob- 
lem can be alleviated by introducing an irregular net with 
smaller spacing near the boundary in such a way that the 
order of magnitude of numerical error is the same throughout 
the flow field. One choice Is to use grids with discontinu- 
ous ly varying resolutions. However, there are two disadvan- 
tages to this method: 1) it 13 necessary to Interpolate 

values of variables or their derivatives at intermediate 
points, and weak numerical instabilities usually arise at 
’•.lie boundary between the large and small grid size, and 2) 
this method cannot give very small grid intervals without 
greatly increasing the computational time. Another 
i>. visibility is to vary the grid intervals continuously, 
avoiding the necessity of intermediate interpolations. This 
can be done by defining a stretched coordinate ; 

♦ * ♦ U) 

such that the grid intervals At are constant while the Ai|t 
varies appropriately, i|> (t) should have the following 
properties : 
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1. ^ should be finite over the whole Interval and 

d? 

Much that there are no singularities in the region of 

interest. If a. becomes Infinite at some point, then the 
dC 

mapping ^ (?) will yield poor resolution near that point, 
which cannot be improved by Increasing the number of points. 



An error analysis based on Taylor series expansion 
suggests the convenience of choosing a function * P n (c), 
where P n is a polynomial of degree greater than one; also, 

P n should be the lowest order polynomial such that good 
resolution is maintained at \p * 0* 

It is also convenient to choose C such that: 

1. Symmetry with respect to C « 0 is restored, and 

2. c coincides with r when the flow is fully developed. 

These requirements strongly suggest an inverse Von 

Mises transformation with fully developed axial velocity j 

replacing axial velocity at any point. For flow of an 
incompressible fluid through a conduit, this yields 

Si - 2 (l-< 2 > C. 

di; 

Solving this differential equation with the boundary condi- 
tion c * 0 at i|> * 0 the following relationship between t|» 

and c is obtained: 

t|/ s ^ r ^ 

2 % 



0 
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It can be shown that this transformation satisfies all 
the above mentioned requirements. The advantages of the c 
transformation are best understood by observing the develop- 
ment of nonunlformly spaced r grids (resulting from uniform- 
ly spaced 5 grids) along the flow. Near the entrance, where 
the gradients are very severe, very small grid sizes are 
maintained close to the wall. As the flow develops, the 
gradients become less severe and r grids become more or less 
uniform, thu3 ensuring errors of the same order of magnitude 
in the computational scheme everywhere In the r direction. 
This also results in a progressive reduction of average 
error in the r direction along the flow. 

A scale factor L is required to map the region of 
interest 0 < r < 1 into 0 < c < 1, in the case of coaxial 
flow of two incompressible fluids. The value of L can be 
shown to be: 


/ o <,u li p ll rdr ^ C U 21 p 21 rdr 
1 (1-r 2 ) rdr 


? 'o 


It Is also convenient to use the same stretch factor 
for the coaxial flow of compressible fluids. 
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APPENDIX B 

ENTRANCE PLOW IN A TUBE 

The results of an analysis on the steady, laminar, 
incompressible, isothermal, Newtonian entrance flow through 
a tube for a uniform entrance velocity were obtained using 
the implicit numerical scheme described in Chapter IV and 
by specifying that: 

a) properties and entrance velocities and temperatures 
of both fluids are the same, and 

b) the heat generation rate is zero. 

The only modification in the numerical scheme is that 
the pressure drop was calculated using a macroscopic momentum 
balance. Results were obtained using a 12 point C grid. 

Even though this is a very coarse grid, the results are in 
good agreement with values obtained by other sources. 

The developing velocity profiles are shown in Figure B-l. 
For comparison purposes, values of dimensionless excess pres- 
sure drop C defined as 

- 2 AP = l6z + C 

and entrance length defined as the dimensionless axial posi- 
tion at which the center line velocity reached 99 % of its 
fully developed value, are listed in Table B-l. 


VEL0CI 






B74 


Table B-l. Values of 0 and Reduced Entrance Lengths 


Source 

C 

Entrance 

Length 

•.’bis .study 

1.24 

0.23 

Campbell and Slattery 1 

1.18 

0.24 

3 

Christiansen and Lemmon 

1.33 

0.24 

ii 

Collins and Schowalter 

1.33 

0.24 

I.mghaar 11 

1.28 

0.23 

Hornbeck^ 

1.28 

0.23 

; rent as ^ 

1.18 

0.23 


The accuracy of the results could further be Improved 
uy increasing the number of grids in the numerical method. 
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